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ABSTRACT 


Precise  navigation  with  high  update  rates  is  essential  for  automatic  landing  of 
an  unmanned  aircraft.  Individual  sensors  currently  available  —  INS,  AHRS,  GPS, 
LORAN,  etc.  —  cannot  meet  both  requirements.  The  most  accurate  navigation  sen 
sor  available  today  is  the  Global  Positioning  System  or  GPS.  However,  GPS  updates 
only  come  once  per  second.  INS,  being  an  on-board  sensor,  is  available  as  often  as 
necessary.  Unfortunately,  it  is  subject  to  the  Schuler  cycle,  biases,  noise  floor,  and 
cross-axis  sensitivity.  In  order  to  design  and  verify  a  precise,  high  update  rate  navi¬ 
gation  system,  a  working  model  of  Differential  GPS  has  been  developed  including  all 
of  the  major  GPS  error  sources  —  clock  differences,  atmospherics.  Selective  Availabil¬ 
ity  and  receiver  noise.  A  standard  INS  system  was  also  modeled,  complete  with  the 
inaccuracies  mentioned.  The  outputs  of  these  two  sensors  —  inertial  acceleration  and 
pseudoranges  —  can  be  optimally  blended  with  a  complementary  Kalman  filter  for 
positioning.  Eventually,  in  the  discrete  case,  the  high  update  rate  and  high  precision 
required  for  autoland  can  be  achieved. 
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I.  INTRODUCTION 


Precise  navigation  has  been  one  of  mankind’s  endeavors  for  thousands  of  years. 
Even  prehistoric  man  had  to  answer  the  question,  “Where  am  I?”  frequently.  Celestial 
navigation,  landmark  navigation,  and  other  forms  of  “natural”  navigation  suffice  for 
some  purposes  even  today.  But  the  human  race’s  non-stop  technological  progress  has 
both  enabled  and  created  a  need  for  more  precise  navigation.  The  advent  of  aviation 
and  space  travel  have  placed  greatly  increased  demands  on  navigation.  Recently 
developed  “smart”  weapons,  as  well,  require  exact  knowledge  of  their  position  in 
order  ..o  fiud  their  targets.  Furthermore,  unmanned  recoverable  aircraft  could  benefit 
greatly  from  improving  navigation  accuracy. 

Today,  work  in  remotely  piloted  vehicljs  is  at  the  forefront  of  technology.  The 
goal  of  the  current  research  and  development  on  the  “Archytas”  (taking  place  here 
at  The  Naval  Postgraduate  School)  is  totally  autonomous  flight  in  all  regimes  — 
automatic  take-off,  transition  from  vertical  take-off  to  horizontal  flight,  automatic 
waypoint  tracking,  and  finally  automatic  landing.  The  task  of  automatically  landing 
an  aircraft,  on  land  or  on  a  pitching  ship  at  sea,  requires  nearly  flawless  guidance, 
navigation,  and  control.  Positioning  errors  which  would  be  considered  minuscule 
under  normal  circumstances  could  easily  be  disastrous  during  an  automatic  landing. 
“Autoland”  retjuires  not  only  exceptional  positioning  accuracy  in  three  dimensions, 
but  also  rapid  updates.  Thus,  the  goal  of  this  undertaking  is  to  provide  the  most 
accurate,  most  rapidly  updated,  three-dimensional  navigation  system  available. 

Many  existing  navigation  systems  provide  either  a  high  update  rate  or  accuracy. 
Probably  the  most  popular  long-range  navigation  system  in  use  is  the  Inertial  Nav¬ 
igation  System  (INS).  This  system  is  entirely  self-contained  on  the  platform  which 


carries  it.  It  senses  aircraft  acceleration  which  it  converts  to  velocity  and  position. 
Inertial  navigation  is  subject  to  numerous  errors,  among  them  the  notorious  Schuler 
cycle,  which  make  the  position  accuracy  degrade  over  time.  This  drift  is  typically 
around  one  nautical  mile  per  hour.  While  the  INS  position  can  be  updated  many 
times  each  second,  its  accuracy  is  insufficient  as  a  stand-alone  navigation  sensor  for 
autoland.  There  are  several  radio  navigation  aids  currently  available  —  GPS,  LO- 
RAN,  OMEGA,  TACAN,  VOR/DME.  They  offer  various  degrees  of  precision  due  to 
their  technical  sophistication  and  overall  age. 

LOR  AN,  short  for  Long  Range  Navigation,  is  a  low  frequency  system  which 
allows  receivers  to  estimate  their  position  by  time-difference-of-arrival.  First,  a  re¬ 
ceiver  calculates  its  position  along  a  hyperbolic,  earth-bound  curve  by  using  a  pair  of 
LORAN  transmitters.  A  second  pair  of  transmitters  allows  the  receiver  to  generate 
another  hyperbola  representing  the  locus  of  possible  positions.  It  is  now  possible 
to  resolve  its  two-dimensional  position  as  the  intersection  of  the  two  hyperbolas. 
LORAN  wavelengths  are  approximately  1.6  miles,  making  truly  precise  positioning 
impossible.  Also,  fixes  are  available  only  ten  to  20  times  per  minute.  This  system  is 
unsuitable  for  autoland  because  of  its  imprecision,  low  update  rate  and  the  fact  its 
fixes  are  only  two-dimensional. [Ref.  1] 

OMEGA  navigation  is  based  on  the  phase-difference-of-arrival  technique. 
Rather  than  using  time  of  arrival,  OMEGA  uses  phase  differences  at  the  receiver 
to  generate  two  hyperbolas,  similar  to  LORAN.  However,  this  system  uses  a  far  lower 
frequency  giving  it  excellent  range.  In  fact,  OMEGA  covers  the  entire  earth  with  only 
eight  transmitters.  This  low  frequency  yields  wavelengths  of  approximately  16  miles. 
Position  accuracy  can  be  expected  to  be  two  to  four  nautical  miles.  LORAN  and 
OMEGA  position  accuracy  can  be  greatly  improved  by  using  the  differential  mode. 
This  mode  allows  a  receiver’s  position  to  be  calculated  relative  to  the  known  position 
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of  another  receiver.  Sub-mile  accuracy  can  be  achieved  with  this  augmented  mode. 
Even  in  the  differential  mode,  the  OMEGA  system  is  unacceptable  for  autoland  as  it 
cannot  provide  the  accuracy  or  the  three  dimensional  fixing  required. [Ref.  1] 

VOR/DME  (VHP  Omnidirectional  Range)  and  TACAN  (Tactical  Air  Naviga¬ 
tion)  systems  both  provide  bearing  and  range  information  to  aircraft  within  their 
line-of-sight.  The  bearing  portion  of  these  systems  is  composed  of  two  signals  — 
a  rotating  “lighthouse”  signal  and  an  omnidirectional  signal.  The  “lighthouse”  por¬ 
tion  of  the  signal  rotates  continuously  at  30  revolutions  per  minute.  Each  time  it 
reaches  magnetic  north,  the  omnidirectional  signal  transmits  one  pulse.  A  receiver’s 
magnetic  bearing  to  the  transmitter  is  a  linear  function  of  the  difference  in  arrival 
time  between  the  omnidirectional  pulse  and  the  “lighthouse  pulse”.  This  method  is 
accurate  to  within  approximately  three  degrees  [Ref.  2,  p.  36].  The  ranging  portion 
of  these  systems  uses  two  way  spherical  ranging.  Equipment  on  board  the  aircraft 
broadcasts  a  pulse  which  the  station  immediately  rebroadcasts  on  another  frequency. 
Once  again,  the  difference  in  time  between  transmission  and  reception  of  this  pulse 
multiplied  by  half  of  the  speed  of  light  yields  slant  range  to  the  station.  Overall  range 
accuracy  of  DME  can  be  expected  to  be  ±0.5  nautical  miles  or  three  parts  in  100, 
whichever  is  greater  [Ref.  Ij.  TACAN  range  accuracy  is  approximately  ±0.1  nautical 
miles.  VOR  uses  VHP  frequencies  while  DME  and  TACAN  use  UHP.  Both  systems 
are  limited  by  line-of-sight.  This  limitation  is  troublesome  for  an  aircraft  needing 
precise  data  at  low  altitude.  While  position  updates  from  these  systems  are  available 
virtually  continuously,  their  accuracies  are  insufficient  for  the  autoland  problem. 

The  remaining  sensor  is  the  Global  Positioning  System  (GPS).  This  is  the  newest 
and  most  accurate  radio  navigation  aid  available.  GPS  provides  positioning  accuracy 
on  the  order  of  100  meters  to  everyone  via  the  Standard  Positioning  Service  (SPS), 
and  16  meters  in  three  dimensions  via  the  Precise  Positioning  Service  (PPS)  to  those 
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authorized  by  the  United  States  Department  of  Defense.  [Ref.  3,  p.  3-1].  GPS 
accuracy  can  be  further  improved  to  just  a  few  meters  by  using  differential  correc¬ 
tions  (see  Chapter  II.  ).  While  its  title  as  most  precise  navigation  aid  available  is 
unchallenged,  its  update  rate  is  currently  once  per  second.  Thus,  it  too  is  insufficient 
as  a  st2Lnd-alone  navigation  sensor  in  the  autoland  phase  of  flight. 

It  seems  that  no  single  sensor  can  provide  highly  precise,  rapidly  updated  posi¬ 
tioning  information.  The  high  update  rate  of  INS  coupled  with  the  superb  precision 
of  differential  GPS  can  solve  the  autoland  navigation  problem.  By  using  the  tech¬ 
niques  of  optimal  estimation,  these  two  sensors  can  be  blended  to  produce  optimal 
estimates  of  position  and  velocity.  The  estimates  will  be  both  precise  and  rapidly 
updated,  thus  meeting  the  established  criteria. 

This  thesis  is  composed  of  4  basic  parts.  The  first  portion  discusses  background 
material  related  to  DGPS/INS  navigation  including  basic  coordinate  systems.  The 
following  section  describes,  in  detail,  the  SiMULiNK  computer  models  of  DGPS  and 
INS.  The  heart  of  this  work,  the  Kalman  filter  which  integrates  the  outputs  of  the 
two  sensors,  is  described  and  verified  in  the  third  section.  The  final  section  offers 
suggestions  for  further  refinements  of  the  concepts  advanced  in  this  thesis. 
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II.  BACKGROUND 


Before  embarking  on  a  study  of  the  intricacies  of  integrated  Global  Positioning 
System  (GPS)  and  Inertial  Navigation  System  (INS)  navigation,  it  is  necessary  to 
be  familiar  with  some  underlying  concepts.  For  example,  a  GPS/INS  navigation 
system  uses  a  minimum  of  three  coordinate  systems.  In  this  chapter,  these  reference 
frames  are  explicitly  defined,  as  well  as  the  transformations  between  them.  Secondly, 
an  overview  of  GPS  is  presented  to  highlight  the  principles  which  are  used  and  the 
errors  which  effect  the  resulting  system.  The  chapter  concludes  with  an  overview  of 
inertial  navigation. 

A.  COORDINATE  SYSTEMS 

GPS-aided  inertial  navigation  involves  the  use  of  three  distinct  coordinate  sys¬ 
tems.  They  are: 

•  earth-centered,  earth-fixed  Cartesian 

•  latitude,  longitude,  altitude  (geodetic) 

•  tangent-plane  Cartesian  (local  geodetic) 

The  earth-centered,  earth-fixed  system  is  independent  of  the  mathematical 
model  of  the  earth’s  surface.  However,  both  the  geodetic  and  the  local  geodetic 
systems  depend  on  the  specification  of  the  earth  model.  The  current  standard  for 
modeling  the  surface  of  the  earth  is  the  WGS-84  ellipsoid.  This  ellipsoid  is  generated 
by  rotating  an  ellipse,  whose  semi-major  axis  is  6378137.0  meters  and  whose  semi¬ 
minor  axis  is  6356752.3  meters,  about  its  minor  axis.  The  resulting  closed  surface  is 
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the  model  of  the  earth’s  surface.  The  true  north  pole  (conventional  terrestrial  pole) 
and  true  south  pole  are  the  endpoints  of  the  minor  axis  of  the  ellipsoid. 

The  remainder  of  this  section  will  define  the  three  coordinate  sy'^tems  in  detail. 

1.  Earth-Centered,  Earth-Fixed  Cartesian  Coordinate  System 

GPS  uses  the  earth-centered,  earth-fixed  Cartesian  system  for  fixing  satel¬ 
lite  positions.  Likewise,  GPS  receivers  calculate  their  navigation  solutions  (t.e.,  po¬ 
sition)  in  the  ECEF  system  before  transforming  to  latitude,  longitude,  and  altitude 
(geodetic).  The  origin  of  this  system,  as  the  name  implies,  is  at  the  center  of  the 
earth.  The  X-axis  goes  through  the  intersection  of  the  equator  (0^^  of  latitude)  and 
the  prime  meridian  (0°  longitude).  The  Y-axis  departs  the  origin  and  passes  through 
the  intersection  of  the  equator  and  90°  E  longitude.  To  complete  the  right-handed 
triad,  the  Z-axis  leaves  the  origin  and  passes  through  the  true  north  pole.  Many 
surveying  oriented  texts  refer  to  this  system  as  the  Conventional  Terrestrial  System 
(CTS).  It  is  depicted  in  Figure  2.1. 


Figure  2.1:  Earth-centered,  Earth-fixed  Coordinate  System 
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2.  Geodetic  Coordinate  System 

The  output  of  navigation  systems  used  on  aircraft  today  is  generally  lati¬ 
tude,  longitude,  and  altitude  —  t.e.  resolved  in  the  geodetic  coordinate  system.  This 
is  the  system  used  for  describing  positions  of  most  earth  bound  objects.  Charts 
developed  for  long  range  land  and  sea  navigation  invariably  use  geodetic  coordinates. 

The  geodetic  coordinate  system  is  somewhat  analogous  to  spherical  coor¬ 
dinates.  The  primary  difference  is  that  the  elevation  angle  or  latitude,  4>,  is  the  angle 
between  the  ellipsoidal  normal  and  the  equatorial  plane.  This  means  that  the  ray 
that  defines  this  angle  does  not  intersect  the  equatorial  plane  at  the  exact  center  of 
the  earth.  Instead,  it  intersects  the  equatorial  plane  at  a  small  radius  outside  of  the 
center  as  shown  in  Figure  2.2.  The  longitude.  A,  is  identical  to  the  spherical  concept 


of  that  angle.  It  is  the  angle  in  the  equatorial  plane  from  0°  latitude  and  0°  longitude 
to  any  given  point.  Finally,  h,  the  geodetic  height  or  altitude,  is  the  distance  along 
the  ellipsoidal  normal  away  from  the  surface  of  the  earth. 
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3.  Tangent  Plane  (Local  Geodetic)  Cartesian  Coordinate  System 

Typically,  pure  inertial  systems  navigate  in  a  so-called  tangent  plane  coor¬ 
dinate  system,  before  outputting  position  in  geodetic  coordinates.  The  tangent  plane 
system  is  defined  by  passing  a  plane  through  any  point  on  the  eairth’s  surface  (see 
Figure  2.3).  The  intersection  of  the  plane  with  the  surface  of  the  earth  becomes  the 


Figure  2.3:  Tangent  Plane  Coordinate  System 


origin  of  the  system.  The  X-axis  points  toward  true  east.  The  Y-axis  points  toward 
true  north.  Lastly,  the  Z-axis  is  perpendicular  to  the  defining  plane  of  the  system, 
away  from  the  center  of  the  earth.  It  is  the  Z  coordinate  of  the  triad  which  defines  a 
point’s  altitude  in  this  system. 

B.  COORDINATE  TRANSFORMATIONS 

In  order  to  use  these  three  coordinate  systems,  one  must  be  able  to  transform  be¬ 
tween  them  freely.  For  example,  satellite  coordinates  enter  the  DGPS  model  (detailed 
in  a  later  chapter)  in  ECEF  coordinates.  However,  the  ranges  to  the  satellites  are 
calculated  in  the  tangent  plane  coordinate  system.  Therefore,  a  transformation  be- 
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tween  ECEF  and  tangent  plane  systems  is  required.  For  other  similar  circumstances, 
transformations  are  necessary  between  all  systems.  With  a  total  of  three  systems,  six 
conversions  are  required  —  one  from  each  system  to  the  two  others,  multiplied  by 
three  systems,  yields  six  transformations.  However,  two  of  these  algorithms  can  be 
assembled  as  combinations  of  the  others.  By  being  able  to  convert  from  geodetic  to 
ECEF,  and  ECEF  to  local  geodetic,  one  can  convert  from  geodetic  to  local  geode¬ 
tic  by  chaining  the  two  conversions  together.  Thus,  the  only  four  transformations 
discussed  are 

•  geodetic  to  ECEF 

•  ECEF  to  geodetic 

•  ECEF  to  tangent  plane 

•  tangent  plane  to  ECEF 

Since  the  vectors  to  be  transformed  are  all  position  (not  free)  vectors,  the  transfor¬ 
mations  are  defined  for  non-free  vectors. 

1.  Geodetic  to  ECEF  Coordinate  Transformation 

In  order  to  compute  the  transformation  from  geodetic  to  ECEF  coordi¬ 
nates,  three  auxiliary  quantities  —  /,  e,  and  N  —  must  first  be  defined.  The  flat¬ 
tening  factor,  /,  represents  the  relative  flatness  of  the  ellipsoid.  A  zero  flattening 
factor  would  mean  an  unflattened  ellipse  (a  sphere),  while  a  unity  value  would  mean 
a  totally  flattened  ellipsoid  (a  circle  in  the  plane  perpendicular  to  the  minor  axis). 
The  mathematical  definition  of  /  is 

/=— ,  (2.1) 

a 

where  a  and  b  are  the  semi-major  and  semi-minor  axes  of  the  ellipsoid,  respectively. 
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Directly  related  to  the  flattening  factor  is  the  eccentricity,  e.  It  is  defined 


by 


=  2/  -  f. 


(2.2) 


The  eccentricity  is  a  variable  similar  to  the  flattening  factor.  It  represents  how  close 
the  ellipsoid  is  to  a  sphere.  It,  too,  is  one  for  a  sphere  and  zero  for  a  completely 
flat  figure.  The  eccentricity,  rather  than  the  flattening  factor,  is  typically  used  in 
coordinate  transformations. 

Lastly,  N  is  the  length  of  the  ellipsoidal  normal  from  the  ellipsoidal  surface 
to  its  intersection  with  the  ECEF  Z  -  axis.  Mathematically,  N  is 


a 

^T— ~?sin^ 


where  <t>  is  the  geodetic  latitude. 

Using  these  quantities,  one  may  define  the  transformation 


(2.3) 


X  =  {N  -^h)  cos<f>  cosA 
y  =  (iV  +  ^)cos^  sinA 

z  =  [iV’(l  —  e*)  +  h]  sin^.  (2.4) 


2.  ECEF  to  Geodetic  Coordinate  Transformation 

The  transformation  from  ECEF  to  geodetic  coordinates  is  clearly  the  in¬ 
verse  of  the  process  presented  in  the  previous  section.  First,  A,  the  longitude,  can  be 
found  by  dividing  first  two  expressions  of  Equations  2.4  yielding 


tanA  =  — . 

X 


(2.5) 


By  examining  the  geometry  in  Figure  2.2,  one  can  determine  the  following  relationship 


tan^  = 


{N  -1-  h)  sm(f> 


(2.6) 
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which  is  a  non-linear  equation  in  <^.  Solving  the  third  of  Equations  2.4  for  {N-\-h)  sin  ^ 
and  substituting  into  Equation  2.6  yields 

e^Nsin(l>, 


tan 


f(l  + 


(2.7) 


which  is  still  an  analytically  unsolvable  equation  in  geodetic  latitude.  To  solve  this 
equation,  one  initially  assumes  that  h  is  zero,  an  excellent  assumption  for  intra- 
atmospheric  flight.  Now,  the  Z  equation  of  Equations  2.4  can  be  simplified  to 


Zh=o  =  N(1  —  e^)sin<^. 


(2.8) 


This  equation  can  be  substituted  into  Equation  2.7  to  give  the  initial  solution  for  <i> 

1 


tan  (f)i  — 


w 


(2.9) 


1  — 

This  initial  solution  for  <f)  can  be  substituted  back  into  the  second  term  of  Equation  2.7 
to  yield  an  updated  <f>.  Iteration  of  this  process  commencing  with  Equation  2.7 
continue  until  the  geodetic  latitude  stops  changing.  Finally,  solving  Equation  2.6 
for  h 

h  =  -  N,  (2.10) 

COS0 

which  completes  the  conversion. 


3.  ECEF  to  Tangent  Plane  (Local  Geodetic)  Coordinate  TVansfor- 
mation 

Anytime  one  uses  a  local  geodetic  or  tangent  plane  coordinate  system,  one 
must  first  specify  the  geodetic  coordinates  —  latitua^  and  longitude  -  of  the  origin. 
Once  specified,  the  origin  must  be  expressed  in  ECEF  coordinates.  Then,  a  vector 
from  this  origin  to  the  point  being  transformed  can  be  formed,  resolved  in  the  ECEF 
system 

{■  Ax  ]  [  X2  -  3:1  ] 

(2.11) 


Ax 

■  X2  -  Xi  ' 

Ay 

= 

y2-y\ 

.  . 

ECEF 

.  ^2  -  . 
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where  the  ‘2’  subscript  denotes  the  point  being  transformed.  The  product  of  two 
rotation  matrices  operates  on  the  difference  vector  defined  in  Equation  2.11  to  yield 
the  local  geodetic  coordinates  of  Pj 


■  I  1 

—  sinA 

cosA 

0  ■ 

Ax 

P2  = 

y 

= 

—  sin(^  cosA 

—  sin<^  sinA 

cos4> 

Ay 

tp 

cos<f>  cosA 

cos<l>  sinA 

sm<f> 

Ac 

where  A  is  the  geodetic  longitude  and  4>  is  the  geodetic  latitude. 


4.  Tangent  Plane  (Local  Geodetic)  to  ECEF  Coordinate  Transfor¬ 
mation 


Unlike  converting  between  geodetic  and  ECEF,  transforming  between  local 
geodetic  and  ECEF  is  invertible.  The  transformation  from  tangent  plane  to  ECEF 
can  be  derived  by  merely  reversing  the  process  developed  in  the  previous  section. 

First,  the  origin  of  the  local  geodetic  system  must  be  converted  from  geode¬ 
tic  to  ECEF  coordinates.  Next,  the  inverse  of  the  rotations  performed  in  Equa¬ 
tion  2.12  must  be  executed  yielding  the  A  vector  from  the  origin  of  the  tangent  plane 
system  to  the  point  being  transformed.  This  vector,  expressed  in  ECEF  coordinates 


'  Ax  ' 

—  sinA  —  sin<^  cosA  cose;!*  cosA 

X 

Ay 

= 

cosA  —  sin<^  sinA  cos(j>  sinA 

y 

(2.13) 

Az 

ECEF 

0  cos4>  sin<^ 

z 

tp 

Clearly,  adding  the  A  vector  to  the  position  of  the  origin  of  the  tangent  plane  system 


(both  now  in  ECEF  coordinates)  completes  the  transformation: 


X 

^orig 

Ax 

y 

— 

yorig 

+ 

Ay 

z 

ECEF 

.  ^orig  , 

ECEF 

Az 

J.  THE  NAVSTAR  GPS 


(2.14) 


The  NAVigation  Satellite  Timing  And  Ranging  Global  Positioning  Sy^stem  is  a 
satellite-based  radio  navigation  system  with  the  capability  to  provide  locating  data 
to  an  unlimited  number  of  users.  The  first  satellite  was  deployed  in  1978  although  the 
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first  receiver  did  not  become  commercially  available  until  1982  [Ref.  4,  section  1.2]. 
This  system  is  the  product  of  experience  gained  from  several  previous  space-based 
navigation  systems  like  TRANSIT  AND  USAF  System  621B.  It  is  comprised  of  three 
segments: 

•  Space  segment 

•  User  segment 

•  Control  segment 

Each  segment  is  discussed  in  the  following  sections. 

1.  Space  Segment 

A  total  of  24  satellites  now  constitute  a  fully  operational  space  segment. 
Twenty-one  of  theses  space  vehicles  (shown  in  Figure  2.4  operate  continuously,  while 
the  remaining  three  act  as  orbiting  spares.  Today,  only  three  Block  I  satellites  remain 
in  orbit.  These  satellites  were  the  first  GPS  satellites  in  space.  The)-  were  launched 


Figure  2.4:  A  NAVSTAR  GPS  Satellite  from  [Ref.  5,  p.  4.01] 
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from  1978  through  1985  [Ref.  4,  section  1.2].  Block  I  space  vehicles  weigh  960  pounds 
and  generate  420  watts  of  electrical  power  [Ref.  2,  p.  27],  The  remaining  satellites  in 
orbit  were  creatively  named  the  Block  II  vehicles.  These  vehicles  are  less  vulnerable 
to  radiation  and  have  more  ntemory.  With  this  increase  in  capability  has  come  an 
increase  in  weight  and  power  requirements.  These  newer  satellites  weigh  2000  pounds 
and  generate  700  watts  of  electrical  pow’er. 

The  space  vehicles  are  inserted  into  orbits  defined  by  the  six  Keplerian 

constants: 


•  semi-major  axis,  a 

•  orbital  eccentricity,  e 

•  orbital  inclination,  i 


•  ascending  node,  fl 

•  argument  of  perigee,  w 


•  time  of  perigee  passage,  T 


The  semi-major  axis  of  a  GPS  satellite  orbit  is  nominally  26,560  kilometers. 
It  is  half  of  the  length  of  the  ellipse  which  defines  the  space  vehicle’s  path.  Second, 
the  eccentricity  (different  from  the  eccentricity  defined  with  respect  to  the  shape  of 
the  WGS-84  ellipsoid)  or  oblateness  of  the  satellites  orbit  is  defined  as  follows; 


rg  ~  rp 

rg+Vp' 


(2.15) 


where  Vg  is  the  apogee  radius  and  Vp  is  the  perigee  radius.  For  GPS  satellites,  the 
eccentricity  cannot  exceed  two  percent.  The  third  of  the  Keplerian  Hements  is  orbital 
inclination.  It  is  the  angle  between  the  plane  defined  by  the  orbit  and  the  equator. 
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For  example,  satellites  in  a  polar  orbit  have  a  90°  orbital  inclination;  those  in  an 
equatorial  orbit  have  0°  orbital  inclination.  Block  I  NAVSTAR  GPS  vehicle  orbits  are 
inclined  at  63°  while  Block  II  satellites  are  inclined  at  55°.  The  ascending  node  is  the 
satellite’s  geodetic  longitude  as  it  passes  through  0°  of  latitude  toward  the  northern 
hemisphere.  GPS  satellites  orbit  in  six  different  planes.  Thus,  there  are  exactly  six 
different  ascending  nodes.  The  last  two  Keplerian  elements  are  the  argument  of  the 
perigee  and  the  t'me  of  perigee  passage.  The  argument  of  the  perigee  is  the  angle 
in  the  orbital  plane  between  the  ascending  node  and  the  closest  point  of  approach 
of  the  satellite  to  earth.  The  time  when  the  vehicle  reaches  this  point  is  the  time  of 
perigee  passage.  The  ranges  of  these  last  two  parameters  span  all  possible  values  for 
GPS  satellites.  That  is,  satellites  reach  their  perigee  at  all  different  times  of  day  and 
different  locations. 

Ideally,  the  six  Keplerian  elements  would  be  sufficient  to  define  any  satel¬ 
lite’s  three-dimensional  position  and  velocity  vectors  for  all  time.  However,  the  ')rbits 
become  perturbed  by  lunar  and  solar  gravity  and  the  earth’s  equatorial  bulge,  &s  well 
as  several  other  less  significant  effects.  Thus,  the  number  of  quantities  required  to 
fully  specify  the  position  and  velocity  of  a  satellite  in  a  real  orbit  is  increased  to  16. 
The  extra  elements  consist  of  time  rates  of  change  of  the  ascending  node  and  inclina¬ 
tion,  i.  e.  n  and  i,  as  well  as  six  other  coefficients  which  account  for  variation  in  the 
earth’s  gravitational  field.  These  16  coefficients  must  become  part  of  the  navigation 
message. 

The  navigation  message  is  the  framework  in  which  the  GPS  satellites  broad¬ 
cast  their  data.  One  “frame”,  the  complete  message,  consists  of  five  subframes  of  300 
bits  each.  Subframe  one  contains  coefficients  used  to  correct  the  satellite’s  clock  to 
exact  GPS  time.  The  16  pieces  of  ephemeris  data  are  broadcast  in  the  second  and 
third  subframes.  Subframe  four  contains  special  messages,  ionospheric  correction  co¬ 
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efficients,  and  coefficients  for  conversion  of  GPS  time  to  Universal  Coordinated  Time. 
Ephemeris  and  health  data  for  the  entire  GPS  constellation  is  transmitted  in  sub¬ 
frame  five.  Because  of  the  volume  of  data  in  subframes  four  and  five,  both  must  be 
subdivided  into  25  pages.  Therefore,  it  takes  25  full  frames,  broadcast  at  the  rate  of 
30  seconds  per  frame,  or  750  seconds  (12.5  minutes)  to  receive  the  entire  navigation 
message.  Critical  navigation  data  -  ephemerides,  and  clock  correction  coefficients  — 
are  updated  every  frame.  Secondary  data  transmitted  in  subframes  four  and  five  is 
provided  primarily  to  assist  the  receiver  in  acquiring  other  satellites.  These  data  are 
not  intended  to  be  precise  so  the  lower  update  rate  of  once  every  12.5  minutes  is 
satisfactory. 

Satellites  use  a  pseudorandom  binary  code  superimposed  on  the  two  carrier 
frequencies  to  communicate.  The  two  frequencies  are  Ll  (1575.42  MHz)  and  L2 
(1227.6  MHz).  The  Ll  frequency  carries  both  the  C/A-  (coarse  acquisition)  and 
P-(precise)  codes,  while  L2  carries  only  the  P-code.  The  less  precise  C/A-code  is 
broadcast  at  a  rate  of  1.023  million  bits  per  second  (Mbps)  and  is  1023  bits  long. 
Therefore,  this  code  repeats  itself  every  millisecond.  The  C/A-code  is  unique  to  each 
satellite  and  does  not  change,  allowing  GPS  receivers  to  quickly  distinguish  between 
space  vehicles,  even  without  access  to  the  P-code.  The  P-code,  being  more  precise, 
is  transmitted  at  10.23  Mbps  with  a  code  length  of  approximately  six  trillion  bits. 
This  code  tcikes  37  weeks  to  repeat.  Since  the  codes  are  reset  every  week  at  midnight 
between  Saturday  and  Sunday,  there  are  sufficient  “code  weeks”  available  in  the  P- 
code  such  that  one  can  be  assigned  to  each  space  vehicle  each  week.  Therefore, 
GPS  receivers  can  easily  distinguish  satellites  from  each  other  with  their  individual, 
weekly-assigned  P-codes.  The  remaining  code  weeks  are  available  for  uplink  from 
the  control  segment  to  the  satellites. 
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2.  User  Segment 

The  many  thousands  of  GPS  receivers  constitute  the  user  segment.  The  re¬ 
ceiver’s  functions  are  to  receive  and  interpret  the  navigation  message  and  to  calculate 
and  output  position.  GPS  receivers  must  determine  the  time  the  navigation  message 
takes  to  travel  the  distance  from  the  satellite  to  the  receiver.  This  is  achieved  by 
autocorrelating  the  pseudorandom  binary  pulse  train  received  from  the  satellite  with 
the  one  in  memory.  A  typical  civilian  GPS  receiver  must  have  the  C/A-codes  for  all 
24  satellites  in  its  memory  (requiring  only  three  kilobytes).  As  the  pseudorandom 
code  is  received,  the  receiver  slews  its  code  until  the  result  of  the  autocorrelation 
function  jumps  to  one.  The  receiver  is  now  “locked-on"  to  that  code.  Multiplying 
the  length  of  time  that  the  receiver  must  slew  its  code  to  achieve  a  unity  correlation 
by  the  speed  of  light  yields  the  “pseudorange” .  This  is  not  the  actual  range  because 
the  offset  of  the  receiver  clock  is  uncertain.  When  the  receiver  can  lock-on  to  four 
satellites  and  thus  measure  pseudoranges  to  each,  is  three-dimensional  position  and 
clock  error  can  be  solved  for  by  inverting  this  set  of  four  equations 

Pi  ~  ^rcvr)^  "I"  {Vsati  Vrcvr)^  (^safi  ~rcvr)^  cAt 

P2  ~  "b  (ysa<2  yrcvr)^  "I”  i^sai2  -^rcvr}^  C^t 

Ps  ~  ^Tcvr)^  "f"  (j/iots  Vtcvt)^  "b  (^»of3  ^rcrr)^  "b 

P4  ~  ^tcvt)^  "b  (j/aat«  J/rcvr)^  "b  ^rcur)^  "b  cAf,  (2.16) 

where  pi  are  pseudoranges,  [xaat,,ysat,,^sat,]  are  the  ECEF  coordinates  of  a  satellite, 
[^rcvr,  ■‘rcvr]  are  the  ECEF  Coordinates  of  the  GPS  receiver,  c  is  the  speed  of 
light,  and  At  is  the  receiver’s  clock  error.  Having  solved  this  set  of  equations,  the 
receiver  now  need  only  transform  the  solution  to  the  geodetic  system  and  display  the 
results. 
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3.  Control  Segment 

The  GPS  Control  segment  is  responsible  for  generating  and  uplinking  clock 
correction  coefficients  and  ephemeris  corrections  for  all  satellites  in  the  constellation. 
Five  control  stations  —  Hawaii,  Ascension  Island,  Diego  Garcia,  Kwajalein,  and  Col¬ 
orado  Springs,  Colorado  -  comprise  this  segment.  These  control  stations  are  essen¬ 
tially  GPS  receivers  capable  of  constantly  tracking  all  satellites  in  view.  Additional 
capabilities  include  highly  accurate  Cesium  clocks  and  recording  facilities.  The  first 
four  of  these  stations  track  all  satellites  in  view  and  record  pseudorange  information 
continuously.  This  data  is  sent  to  the  Master  Control  Station  at  Colorado  Springs 
where  it  is  processed.  When  all  four  monitor  stations  have  locked-on  to  a  single 
space  vehicle  simultaneously,  the  exact  inverse  of  the  standard  navigation  problem 
mentioned  in  the  previous  section  exists.  Since  the  exact  locations  of  the  monitor 
stations  are  known  and  their  clocks  are  extremely  accurate,  the  four  unknowns  be¬ 
come  the  three  coordinates  of  the  space  vehicle  position  and  its  clock  offset.  The 
Master  Control  Station  calculates  these  quantities  and  from  them,  derives  the  nec¬ 
essary  ephemeris  and  clock  corrections.  This  information  is  uplinked  to  each  space 
vehicle  at  least  daily. 

4.  Differential  GPS 

Although  GPS  alone  provides  highly  accurate  positioning,  it  can  be  made 
still  more  accurate  by  augmenting  it  with  a  differential  station.  A  differential  station 
is  merely  another  GPS  receiver  whose  exact  location  is  known.  When  this  second 
receiver  is  near  the  first  receiver,  both  are  subject  to  nearly  the  same  errors,  i.  e.  the 
S2ime  local  atmospheric  properties,  nearly  identical  elevation  angles  and  propagation 
paths  to  any  given  GPS  satellite,  the  same  clock  errors  and  the  same  ephemeris 
errors  for  each  satellite.  By  employing  the  Pythagorean  theorem  on  its  position  and 
the  satellite  position  broadcast  in  the  navigation  message,  the  differential  station 
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can  calculate  the  exact  range  to  that  satellite.  Meanwhile,  it  can  also  calculate 
pseu  ange  in  the  standard  way  (see  Section  2.).  By  comparing  these  two  values 
for  each  satellite  in  view,  the  differential  station  can  evaluate  the  pseudorange  error 
to  each  satellite.  These  values  can  be  broadcast  periodically  to  be  used  by  receivers 
in  the  local  area  to  improve  accuracy.  By  using  the  differential  corrections,  GPS 
accuracy,  even  for  single  frequency  C/A-code  users,  can  be  improved  to  one  to  seven 
meters  rms  [Ref.  2,  p.  64]. 

D.  GPS  ERROR  SOURCES 

Despite  its  exceptional  accuracy,  GPS  is  subject  to  numerous  error  sources. 
Clearly,  the  major  error  sources  must  be  included  in  the  DGPS  model.  Error  sources 
are: 

•  atmospheric  delays 

•  Selective  Availability 

•  clock  differences 

•  ephemeris  error 

•  multipath 

•  receiver  noise 

•  Dilution  of  Precision  (DOP) 

Each  of  these  error  sources  is  ;K,scussed  in  detail  in  the  following  sections. 

1.  Atmospheric  Delays 
a.  Ionospheric  Delays 

The  ionosphere  is  a  layer  of  charged  particles  between  100  and  1000 
kilometers  above  the  earth’s  surface.  These  particles  interact  with  the  transmitted 
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GPS  signal  and  slow  it,  increasing  pseudoranges.  The  equation  describing  this  delay 
is 

40.3 

At  =  —  TEC,  (2.17) 

where  At  is  the  delay  in  seconds,  c  is  the  speed  of  light  (3  x  10*  m/s),  /  is  the 
system  frequency  (1575.42  MHz  for  LI),  and  TEC  is  the  Total  Electron  Content 
(electrons/m^)  along  the  signal's  path.  The  TEC  is  strongly  effected  by  the  solar 
cycle,  season,  time  of  day,  and  latitude.  TEC  maxima  occur; 

•  daily:  1400  local 

•  setisonally:  spring  equinox 

•  solar  cycle;  every  11  years  (next  2001-2002) 

•  geographic:  20°  magnetic  latitude 

It  varies  ±  25%  rms  at  all  latitudes  during  daylight  [Ref.  4,  section  2.5,  p.  24].  The 
pseudorange  error  due  to  ionospheric  effects  can  be  as  great  as  40  meters  [Ref.  4, 
section  2.5,  p.  13]. 

The  algorithm  described  in  [Ref.  6]  which  removes  55-60%i  of  the 
ionospheric  delay  is  based  on  Figure  2.5,  which  shows  the  typical  diurnal  variation  of 
the  ionospheric  delay.  The  “ACTUAL  DATA”  curve  shown  in  Figure  2.5  is  modeled 
with  the  “COSINE  MODEL”  curve,  a  half-cosine.  The  equation  of  the  model  curve 
is 

At  =  DC  +  A  (2.18) 

where  DC,  A,  Tp,  and  P  (constant  offset,  amplitude,  phase,  and  period,  respectively) 
describe  the  diurnal  variation  of  the  ionospheric  delay.  DC  and  Tp  are  assumed 
constant  at  five  nanoseconds  and  1400  local  time,  respectively.  Amplitude  and  period 
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Figure  2.5:  Diurnal  Ionospheric  Delay  from  [Ref.  6,  p.208] 


are  each  modeled  as  four  term  power  series  as  follows: 

^ 

n=0 

P  =  ilMl,  (2.19) 

n=0 

where  the  On’s  and  /?„’s  are  constants  which  are  broadcast  in  the  GPS  navigation 
message,  chosen  based  on  the  day  of  the  year  and  average  solar  flux  over  the  past  five 
days,  and  <f>m  is  the  geomagnetic  latitude  of  the  ionospheric  subpoint.  The  ionospheric 
subpoint  is  the  intersection  of  the  line  between  the  space  vehicle  and  the  receiver  with 
the  surface  at  the  mean  height  of  the  ionosphere. 

Next,  one  must  first  find  the  subtended  earth  angle  between  the  user 
and  the  satellite,  EA  (degrees) 

445 

EA  =  -  4,  (2.20) 

e/  +  20  ’  ^ 

where  el  is  the  elevation  angle  of  the  satellite  with  respect  to  the  user  in  degrees. 
Knowing  EA,  the  geodetic  location  of  the  ionospheric  subpoint  can  be  approximated 

by 


<i>I  =  </'rci;r  +  ^^AcOSa^ 


^  EAsinaz 

^TCVT  H  » 

COS  <1>I 


(2.21) 


A/  = 


where  <f>  and  A  denote  geodetic  latitudes  and  longitudes,  respectively  and  az  is  the 
azimuth  of  the  satellite  with  respect  to  the  receiver.  Now  the  geodetic  latitude  can 
be  converted  to  geomagnetic  latitude  (the  required  quantity)  with  the  following  ap¬ 
proximation: 

=  <^;-|-11.6cos(A;-  291),  (2.22) 

where  all  angles  are  in  degrees. 

The  dimensionless  scale  factor  (SF)  which  scales  the  entire  delay  is 
=  (2-23) 

The  final  expression  for  At  is  a  three  term  Taylor  series  expansion  of 
Equation  2.18  is 

At 
where 

X 

Due  to  the  25%  rms  variation,  the  error  is  modeled  with  a  25%  stan¬ 
dard  deviation.  Therefore,  the  random  part  of  the  error  remains  within  ±25%^  of 
nominal  68%  of  the  time. 

b.  Tropospheric  Delays 

The  lower  section  of  the  atmosphere  also  causes  signal  propagation 
delays.  Typical  tropospheric  delay  is  approximately  two  meters  for  90°  satellite  el¬ 
evation  (directly  overhead)  up  to  28  meters  at  a  five  degree  elevation  angle  [Ref.  7, 
p.  218].  In  this  application,  the  atmosphere  can  be  modeled  as  being  composed  of 


■( 


SF-[DC-¥A{1-^  +  0  for 
SF-[DC]  ■  for 


l<! 


2t(<  -  Tp) 


(2.24) 
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“wet  air”  and  “dry  air”.  Dry  air  is  responsible  for  90%  of  the  total  tropospheric  de¬ 
lay,  whereas,  wet  air  is  responsible  for  only  ten  percent.  While  the  moisture  content 
in  the  troposphere  is  virtually  impossible  to  model  accurately,  this  inaccuracy  has 
minimal  impact.  Numerous  models  which  calculate  the  tropospheric  delay  have  been 
developed.  Black  developed  the  following  model  in  [Ref.  8].  Let 


ajid  As  is  the  wet  or  dry  delay  in  meters,  r,  is  the  distance  from  the  center  of  the  earth 
to  the  station,  P  is  the  surface  pressure  in  atmospheres,  E  is  the  satellite  elevation 
angle,  and  T  is  the  surface  temperature  in  degrees  Celsius.  It  should  be  noted  that  4 
is  an  empirical  constant.  The  value  of  0.85  is  only  valid  for  elevation  angles  above  five 
degree  (GPS  receivers  typically  ignore  satellites  at  lesser  elevation  angles).  Similarly, 
kyj  is  an  empirical  constant  which  varies  based  on  latitude  and  season.  The  value  0.20 
corresponds  to  the  value  for  spring  or  fall  in  mid-latitudes. 
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This  model  has  been  shown  to  be  virtucdly  exact  at  elevation  eingles 
greater  than  40°,  with  its  worst  error  of  about  0.045  m  occurring  between  five  and 
ten  degrees  of  elevation. 

This  tropospheric  model  is  assumed  to  vary  15%  from  the  nominal 
value.  Therefore,  it  is  modeled  with  a  7.5%  nominal  standard  deviation.  This  main¬ 
tains  the  random  part  of  the  error  within  15%  of  the  model  value  95%  of  the  time. 

2.  Selective  Availability 

Selective  Availability  is  a  method  that  the  Department  of  Defense  can  use 
to  intentionally  degrade  the  accuracy  of  pseudorange  measurements  Typically,  this  is 
accomplished  by  dithering  the  space  vehicle  clock  signal.  Dithering  the  clock  involves 
encoding  the  binary  time  signal  the  space  vehicle  broadccists.  The  decryption  process 
is  classified  and  available  only  to  DOD  authorized  users. 

The  use  of  SA  essentially  results  in  the  satellites’  “lying”  to  the  receiver 
about  their  position.  Clearly,  this  adversely  impacts  precision.  Currently,  SA  is  in  op¬ 
eration  on  all  Block  II  space  vehicles  which  comprise  the  majority  of  the  constellation. 
The  DOD’s  stated  goal  for  the  positioning  accuracy  under  SA  is  100  meters  (twice 
rms)  for  a  two-dimensional  fix  [Ref.  9].  According  to  [Ref.  10,  p.  420],  the  selective 
availability  error  can  be  modeled  as  a  zero-mean,  five  meter  standard  deviation  low 
frequency  noise.  The  suggested  cutoff  frequency  for  SA  noise  is  1/180  Hz. 

3.  Clock  Differences 

The  clock  model  used  in  this  treatment  of  DGPS  is  a  two  state  model 
shown  in  Figure  2.6.  It  is  reasonable  to  expect  the  clock  to  have  both  a  bias  and 
drift.  From  daily  exposure  to  clocks,  the  average  person  realizes  that  most  clocks 
are  slightly  offset  from  correct  time  (bias),  and  that  their  accuracy  tends  to  degrade 
with  time  (drift).  From  an  engineering  standpoint,  these  two  phenomena  can  be  best 
modeled  with  zero  mean,  white,  Gaussian  noise.  Rewriting  the  model  in  state  space 
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Figure  2.6:  Receiver  Clock  Model 


form  for  further  analysis  yields 


where 


r/  T\  [  <51  0 

>  =  [  0  5a  ■ 


(2.26) 


The  covariance  of  the  clock  error  state  (xj)  can  be  found  by  solving  the  Lyapunov 
equation  which  can  be  found  numerous  control  textbooks,  one  of  which  is  [Ref.  11,  p. 
104].  This  equation  must  be  solved  over  a  finite  time  interval  (At)  because  the  two 
state  clock  model  is  unstable.  This  interval  would  normally  be  the  sampling  time,  if 


the  model  were  discrete.  The  result  is 


E{xl)  =  SiAt  + 


(2.27) 


By  taking  the  square  root  of  the  variance  and  dividing  by  At,  one  finds  the  more 


standard  clock  parameter  the  Allan  variance,  Va 

I  Si  52  At 

''•'  =  VS  +  — • 


(2.28) 


A  representative  plot  of  the  two  state  model  Allan  variance  as  a  function  of  averaging 
time  (i.  e.  At)  is  shown  in  Figure  2.7. 
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Figure  2.7:  Ideal  Allan  variance 


Real  clocks  behave  somewhat  differently  from  this  simple  model.  A  typical 


Allan  variance  plot  for  a  real  clock  is  shown  in  Figure  2.8. 


Figure  2.8:  Real  Allan  variance 
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The  flat  portion  of  the  curve  is  called  the  flicker  floor.  It  is  the  result  of 
a  non-linear  effect  which  cannot  be  modeled  by  the  two  state  model.  This  causes  a 
significant  discrepancy  between  this  simple  model  and  the  real  world. 

In  order  for  the  model  to  better  represent  reality,  it  must  be  carefullj 
cr2dted  to  fit  the  actual  plot  as  much  as  possible.  By  carefully  choosing  the  v<ilues  of 
Si  and  52,  it  is  possible  to  make  the  actual  and  model  curves  fairly  close.  The  key 
Allan  variance  parameters  between  0.1  and  ten  seconds  of  averaging  time  are  ho,  h-i, 
and  /i_2  [Ref.  12],  Values  of  these  three  parameters  for  three  common  GPS  timing 
standards  are  shown  in  the  Table  2.1  from  [Ref.  10,  p.  428). 

TABLE  2.1:  ALLAN  VARL4NCE  PARAMETERS  FOR  THREE  COM¬ 
MON  TIMING  STANDARDS 


timing  standard 

ho 

h-i 

■M 

crystal 

[Qsiii 

ovenized  crystal 

lOBOSII 

0 

X 

Rubidium 

IBBlHi 

BESmi 

QPQggil 

Brown  [Ref.  10]  finds  that  one  must  choose  where  the  2-state  model  is 
accurate.  Since  normal  averaging  times  are  in  the  0.1  to  ten  second  interval  already 
mentioned,  this  is  the  region  where  the  model  is  made  accurate.  Maximizing  the 
accuracy  of  the  model  in  that  region  dictates  the  following  values  for  the  noise  spectral 
densities: 

ho 

rs,/  — — 

2 

~  2n^h.2  s"^  (2.29) 

In  order  to  remain  conservative  in  this  generic  model  of  DGPS,  the  least 
accurate  clock  —  the  crystal  clock  —  is  used.  The  values  for  Si  and  S2  for  this  clock 
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are 


51  =  4  X  10-*^ 

52  =  1.58  X  10-‘®  s-^  (2.30) 

These  values  are  be  used  in  the  DGPS  error  model. 

4.  Ephemeris  Error 

In  converting  the  pseudoranges  of  at  least  four  satellites  (six  in  this  model) 
to  a  three-dimensional  position  aind  clock  error,  one  must  solve  a  series  of  non-linear, 
coupled  algebraic  equations.  In  these  equations,  the  positions  of  the  satellites  are 
critical.  The  only  way  a  GPS  receiver  or  navigation  filter  knows  the  space  vehicle 
positions  is  the  navigation  message.  If  broadcast  satellite  positions  are  incorrect, 
the  accuracy  of  the  resulting  receiver  position  suffers.  The  control  segment  of  the 
GPS  system  maintains  positions  on  the  entire  constellation  quite  accurately.However, 
it  would  be  folly  to  expect  the  satellites  to  broadcast  inerrantly  accurate  positions. 
Typical  ephemeris  inaccuracies  according  to  [Ref.  13]  are  shown  in  Figure  2.9. 

These  errors  Me  resolved  in  a  coordinate  system  local  to  each  space  vehi¬ 
cle.  The  three  mutually  perpendicular  directions  are  radial,  along  track,  and  cross 
track.  Because  the  DGPS  model  developed  in  this  thesis  does  not  account  for  satel¬ 
lite  motion,  along  track  and  cross  track  directions  cannot  be  resolved.  To  remain 
conservative,  these  two  errors  are  combined  into  a  circular  error  in  the  plane  they 
define.  This  error  is  modeled  as  zero  mean,  white,  Gaussian  noise  with  a  standard 
deviation  of  3.5  meters.  Likewise,  the  radial  error  is  modeled  as  zero  mean,  white 
Gaussian  noise  with  a  0.5  meter  standard  deviation.  Both  of  these  errors  make  the 
stated  accuracy  the  two  rms  point.  In  other  words,  the  value  stays  within  the  stated 
accuracy  95%  of  the  time. 

Since  the  space  vehicle  location  enters  the  model  in  earth  centered,  earth 
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fixed  Cartesian  coordinates  while  the  error  is  added  in  geodetic  coordinates,  a  trans¬ 
formation  is  between  the  two  coordinate  systems  must  be  performed  (see  Subsection  2. 
of  Section  B.,  Chapter  II.  ).  Having  converted  the  satelTite  positions  to  geodetic  co¬ 
ordinates,  the  random  errors  can  be  added.  However,  since  the  geodetic  coordinates 
contain  amgles,  the  random  errors  in  position  must  be  converted  to  equivalent  an¬ 
gles  in  latitude  and  longitude  by  dividing  by  the  radius  of  the  space  vehicles’  orbits, 
~26,560,000  meters.  The  position  and  the  error  are  now  in  compatible  coordinates 
and  can  thus  be  summed. 

5.  Multipath 

The  signal  radiated  by  a  satellite  is  not  required  to  take  a  direct  path  to 
a  receiver.  If  the  signal  encounters  an  electromagnetically  reflective  object,  it  may 
bounce  off  that  object  juid  still  find  its  way  to  the  receiver.  This  signal  has  now 
traveled  a  greater  distance  than  the  straight  line  joining  the  space  vehicle  and  the 
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receiver.  Because  the  receiver  assumes  the  direct  path  is  used,  this  multipath  phe¬ 
nomenon  can  introduce  pseudorange  errors.  Due  to  the  satellite/receiver  geometry, 
multipath  is  far  more  likely  at  low  elevation  angles.  The  GPS  system  has  several 
attributes  that  minimize  multipath  effects  [Ref.  4]: 

•  The  L  band  frequency  (1227.6  MHz)  tends  to  undergo  diffuse  rather  than  spec¬ 
ular  reflection. 

•  The  receiver  antennas  tend  to  reject  multipath  signals. 

•  The  navigation  message  is  broadcast  with  circular  polarization.  Circularly  po¬ 
larized  signals  undergo  reversal  upon  reflection. 

•  GPS  receivers  generally  use  mask  angles  (the  elevation  angle  below  which  the 
satellite  is  ignored)  of  at  least  five  degrees. 

All  of  these  factors  tend  to  attenuate  the  strength  of  any  reflected  signal  making  the 
multipath  effect  insignificant.  Therefore,  it  is  not  modeled. 

6.  Receiver  Noise 

All  receivers  corrupt  the  signals  they  receive.  GPS  receivers  are  no  excep¬ 
tion.  Inaccuracies  resulting  from  quantization  error,  loop  tracking  errors  and  other 
hardware  inadequacies  corrupt  the  pseudorzinge  accuracy.  According  to  [Ref.  4],  rep¬ 
resentative  GPS  receiver  noise  has  a  standard  deviation  of  7.5  m.  Receiver  noise  is 
thus  modeled  as  a  zero  mean,  white  Gaussian  noise  with  7.5  m  standard  deviation.  .4s 
the  differential  station  and  the  aircraft  receivers  are  eissumed  identical,  both  aircraft 
pseudoranges  and  differential  station  pseudoranges  are  contaminated  with  this  noise. 
The  noise  in  the  two  receivers  is  assumed  to  be  independent. 
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7.  Dilution  of  Precision 

All  of  the  errors  discussed  to  this  point  directly  effect  pseudorange  accu¬ 
racy.  The  various  delays  and  noise  sources  can  cause  the  receiver’s  evaluation  of  range 
to  the  space  vechicle’s  to  be  inaccurate.  However,  the  pseudoranges  themselves  are 
irrelevant.  The  science  of  navigation  is  concerned  with  positioning.  Dilution  of  Pre¬ 
cision  is  the  effect  that  links  pseudorange  accuracy  to  position  accuracy.  DOP  can 
be  further  classified  into  several  different  types: 

•  VDOP  —  Vertical  DOP  (z) 

•  HDOP  —  Horizontal  DOP  (x,y) 

•  PDOP  —  Position  DOP  (x,y,z) 

•  TDOP  —  Time  DOP  (t) 

•  GDOP  —  Geometric  DOP  (x,y,z,t) 

The  equation  which  relates  DOP  and  pseudorange  error  is 

<7p  =  DOP  •  <7o,  (2.31) 

where  <Tp  is  the  standard  deviation  of  the  position  error  and  Oq  is  the  standard  devia¬ 
tion  of  the  pseudorange  error,  often  called  the  User  Equivalent  Range  Error  (UERE). 

Dilution  of  Precision  is  a  function  of  satellite  to  receiver  geometry.  As 
Figure  2.10  shows  for  a  four  satellite  constellation,  GDOP  is  minimized  wdth  the 
space  vehicles  spread  out  as  much  as  possible.  In  fact,  the  volume  of  the  tetrahedron 
formed  by  the  unit  vectors  from  the  receiver  to  each  satellite  is  an  empirical  measure 
of  DOP.  PDOP  is  inversely  proportional  to  the  volume  of  the  tetrahedron.  If,  for 
example,  all  of  the  space  vehicles  a  receiver  was  using  for  positioning  were  in  the 
same  plane,  PDOP  would  approach  infinity  (the  volume  of  the  tetrahedron  would  be 


31 


POOR  GDOP  GOOD  GDOP 

satellites  bunched  (ideal  case) 

together  *  one  satellite  overhead 

•  3  on  horii;on, 

1211°  apart  in  aziniiith 

Figure  2.10:  Dilution  of  Precision  from  [Ref.  5,  p.  4.22] 

zero).  Likewise,  minimum  PDOP  is  achieved  with  the  geometry  shown  at  the  right 
of  Figure  2.10. 

The  goal  for  the  design  of  the  entire  GPS  constellation  of  satellites  is  a 
PDOP  no  greater  than  six  everywhere  on  the  earth.  With  the  entire  set  of  24  space 
vehicles  now  in  orbit,  users  can  expect  PDOP  values  under  three  [Ref.  14]. 

The  entire  concept  of  GPS  navigation  heis  now  been  thoroughly  discussed. 
All  of  the  information  put  forth  in  this  discussion  will  be  used  in  developing  the  DGPS 
computer  model  in  the  next  chapter.  To  complete  the  sensor  discussion,  a  description 
of  INS  follows. 


E.  INERTIAL  NAVIGATION 


Inertial  navigation  has  long  been  the  steindard  for  self-contained,  long-range 
aircraft  navigation.  An  Inertial  Navigation  System  (INS)  senses  aircraft  thrust  accel¬ 
eration,  angular  rates  and  spatial  orientation  resolved  in  an  orthogonal  system  and 


computes  the  inertial  acceleration.  This  acceleration  can  now  be  integrated  —  once 
to  find  velocity,  and  twice  to  find  position. 

The  primary  component  in  <uiy  inertial  navigation  system  is  the  Inertial  Measur¬ 
ing  Unit  (IMU).  The  IMU  is  composed  of  three  accelerometers,  three  rate  gyros,  and 
two  inclinometers.  The  accelerometers  measure  thrust  acceleration.  Thrust  accelera¬ 
tion  is  composed  of  lineiu:,  centripetal,  and  gravitational  effects.  Einstein’s  Principle 
states  that  it  is  impossible  for  a  sensor  to  distinguish  between  the  effects  of  gravity 
and  acceleration.  Thus,  the  thrust  acceleration  that  it  provides  is 

=^v  +  X  +\  (2.32) 

where  ®a  is  the  thrust  acceleration,  is  the  linear  acceleration,  is  the  angular 
velocity,  is  the  velocity,  and  is  gravity.  All  quantities  axe  in  the  body  frame. 
This  principle  necessitates  that  the  computation  portion  of  the  INS  compute  and 
remove  local  gravity  from  the  measured  accelerations. 

Ra.te  gyros  sense  angular  velocities.  These  angular  velocities  can  be  resolved 
in  either  the  body  or  inertial  frame,  depending  on  the  type  of  IMU  implementation. 
The  two  inclinometers  sense  aircraft  inertial  orientation,  z.e.,  Euler  angles. 

There  are  two  conceptual  methods  for  implementing  inertial  navigation: 

•  Gimbaled  IMU 

•  Strapdown  IMU 

Brief  discussions  of  each  method  follow. 

1.  Gimbaled  IMU 

A  gimbaled  IMU  rotates  about  its  four  gimbals  during  operation  (see  Fig¬ 
ure  2.11.  Aircraft  IMU’s  must  have  four  gimbals  to  prevent  gimbal  lock,  while  earth- 
bound  IMU’s  require  only  three  gimbals  [Ref.  15,  p.  192].  A  controller  maintains 
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the  IMU  in  a  constant  inertial  orientation  toward  true  North  and  lying  in  the  locally 
horizontal  plane.  By  maintaining  this  orientation,  the  gimbaled  IMU  measures  iner¬ 
tial  quantities  directly.  This  data  can  be  integrated  without  transformation  to  yield 
inertial  velocity,  position  and  Euler  angles.  From  a  navigation  standpoint,  this  con- 
hguration  seems  ideal.  However,  gimbaled  systems  are  large  and  heavy,  making  them 
impractical  for  small  aircraft.  Also,  kinematic  quantities  resolved  in  the  aircraft -fixed 
coordinate  system,  necessary  for  stability  augmentation  ajid  control,  are  not  directly 
available.  Instead,  they  must  be  computed  through  a  series  of  Euler  rotations.  The 
computations  required  take  time  and  thus  introduce  delays  into  an  often  time  critical 
control  problem.  This  fact  makes  the  gimbaled  system  less  than  desirable  for  control. 

2.  Strapdown  IMU 

Strapdown  IMU  is  conceptually  the  reverse  of  the  gimbaled  system  men¬ 
tioned  above.  Rather  than  maintaining  a  constant  inerti€d  orientation,  a  strapdown 
system  is  “strapped  down”  to  the  aircradt,  thus  maintaining  a  constant  orientation 
in  the  aircraft-fixed  coordinate  system.  Therefore,  the  output  of  the  IMU  is  resolved 
in  the  local  coordinate  system.  Kinematic  quantities  are  immediately  available  for 
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the  control  system.  The  extra  computational  burden  now  rests  on  the  navigation 
computer  which  must  transform  the  accelerations  and  angular  velocities  sensed  in 
the  local  coordinate  system  to  the  inertial  system.  Currently,  most  inerticil  systems 
are  of  the  strapdown  variety.  With  the  advent  of  high  speed,  low  cost,  lightweight 
computing  power,  the  required  calculations  in  transforming  from  the  aircraft -fixed 
to  the  inertial  coordinate  system  are  no  obstacle  to  navigation.  This  is  the  system 
which  is  modeled  in  this  thesis. 


F.  INS  COMPUTATIONS 


In  the  strapdown  configuration,  the  IMU  measures  angular  velocities  and  thrust 
acceleration  in  the  body  frame,  as  well  as  Euler  angles.  However,  the  INS  must  provide 
position  and  orientation,  both  in  the  inertial  frame.  In  order  to  compute  position  and 
orientation  in  the  inertial,  tangent  plane  coordinate  system,  inertial  accelerations  and 
Euler  rates  must  be  calculated. 


Before  converting  the  ine’^tial  acceleration  from  the  body  frame  to  the  inertial 
frame,  it  is  necessary  to  compute  the  inertial  orientation.  Computing  the  Euler  rates 
is  a  tricky  endeavor.  The  Euler  rates  are  simply  related  to  the  body  angular  rates  by 
Poisson’s  equation 
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where  p,  q,  and  r  are  the  components  of  ®u;,  the  body’s  angular  velocity.  Obvious!}-, 
this  formula  requires  the  exact  Euler  angles.  However,  one  can  only  measure  these 
angles  directly  at  very  low  frequency.  At  frequencies  greater  than  a  fraction  of  a  Hertz, 
one  must  integrate  the  Euler  rates  found  from  Equation  2.33  to  find  the  angles.  This 
presents  the  seeming  paradox  of  needing  to  know  the  Euler  angles  in  order  to  find 
the  Euler  rates  which  must  be  integrated  to  find  the  Euler  angles.  This  process 
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must  be  implemented  recursively.  That  is,  the  results  of  the  integration  to  find  the 
angles  must  be  fed  back  into  Equation  2.33.  Furthermore,  a  complementary  Kalman 
filter  is  necessary  to  provide  an  optimal  estimate  of  the  Euler  angles,  trusting  the 
inclinometers  at  the  low  frequencies  and  the  rate  gyros  at  high  frequencies.  This 
process  will  be  discussed  in  detail  in  Chapter  IV. 

Now  having  inertial  acceleration  in  expressed  in  the  body  frfime  and  the  orien¬ 
tation  of  the  body  (Euler  angles),  the  transformation  from  body  to  inertial  can  be 
executed.  This  coordinate  transformation  is  defined  by 

+  (2.34) 

where  $,  0,  and  are  the  roll,  pitch  and  yaw  Euler  angles,  respectively,  ^'P  is  the 
aircraft’s  acceleration  in  inertial  coordinates,  is  thrust  acceleration  in  the  aircraft- 
fixed  coordinate  system,  is  gravity  in  inertial  coordinates  and  gR  is  the  transfor¬ 
mation  matrix  from  the  body-fixed  coordinates  to  inertial  tangent  plane  coordinates 
as  follows 

cos  cos  0  cos  sin  0  sin  $  —  sin  $  cos  $  cos  'I'  sin  0  cos  $  -f  sin  $  sin  $ 

—  sin  ♦  cos  0  —  sin  0  sin  $  sin  ^  -I-  cos  cos  $  -  sin  0  cos  $  sin  -  cos  4'  sin  $ 

sin  0  —  cos  0  sin  $  —  cos  0  cos  $ 

It  is  which  can  be  integrated  to  provide  velocity  and  position  in  the  inertial, 

tangent  plane  frame. 

G.  ms  ERROR  SOURCES 

The  Inertial  Measuring  Unit  is  subject  to  a  few  main  error  sources.  These  are: 

•  bias 

•  cross-axis  sensitivity 

•  noise  floor 
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All  of  these  errors  apply  to  both  accelerometers  and  rate  gyros.  Brief  descriptions  of 
these  problems  follow. 

1 .  Biases 

A  bias  in  an  accelerometer  or  a  rate  g>’ro  is  defined  as  a  constant  offset  of 
the  output  of  the  device  from  the  true  value.  In  other  words,  an  accelerometer  might 
constantly  read  0.01  m/s^  while  the  device  is  not  accelerating.  This  0.01  m/s^  would 
be  referred  to  as  a  bias.  This  error  is  modeled  with  a  series  of  small  step  functions, 
one  for  each  accelerometer  and  rate  gyro. 

2.  Cross-Axis  Sensitivity 

Cross-axis  errors  are  caused  by  misalignment  of  the  IMU  with  the  aircraft 
coordinate  axes.  Ideally,  the  components  of  the  IMU  —  the  accelerometers  and  the 
rate  gyros  —  would  each  be  perfectly  aligned  with  the  three  aixes  of  the  aircraft-fixed 
coordinate  system.  Unfortunately,  even  the  highest  fidelity  inertial  sensors  are  never 
precisely  coincident  with  the  appropriate  axes.  Misalignment  of  both  types  of  devices 
causes  errors.  The  cross-axis  errors  are  modeled  with  the  following  equation: 
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where  a^a  is  the  cross-axis  error,  ex,  ey,  and  ez  are  the  cross-axis  error  terms  and 
a  is  the  actual  acceleration  resolved  in  the  aircraft  coordinate  system.  The  e’s  from 
the  previous  equation  are  determined  by  the  amount  of  angular  offset  of  each  sensor 
from  the  correct  position. 

3.  Noise  Floor 

All  sensor  devices  corrupt  the  quantities  they  measure.  Various  factors 
including  thermal  noise  can  cause  a  constant  output  of  white  noise,  regardless  of  the 
actual  acceleration  (or  angular  velocity).  This  noise  floor  makes  accelerations  below  it 
not  measurable.  That  is,  the  output  of  the  sensor  still  includes  the  actual  quantities. 
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but  it  is  “invisible”  because  the  noise  floor  obscures  it.  This  process  is  modeled  by  in¬ 
troducing  a  threshold  to  the  actual  acceleration  in  the  aircraft  coordinate  system  and 
adding  white  noise.  Taking  these  steps  result  in  the  accelerometer  always  reporting 
noisy,  zero-mean  acceleration  unless  the  actual  value  is  above  the  threshold.  VV’hen 
the  actual  value  exceeds  the  threshold,  that  value  is  added  to  the  noise  floor  value 
(and  the  bias)  to  create  the  output  of  the  accelerometer.  This  process  also  applies  to 
the  rate  gyros. 

Now  both  of  the  sensors  —  GPS  and  INS  —  have  been  completely  inves¬ 
tigated.  In  the  following  chapter,  computer  models  of  both  sensors  are  constructed 
using  the  basic  principles  of  operation  of  each  sensors  and  all  of  the  error  sources. 
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III.  SiMULiNK  MODEL 


A.  DGPS  SiMULiNK  MODEL 

The  discussion  of  the  basic  principles  of  GPS  navigation  and  the  various  phe¬ 
nomena  effecting  its  performance  can  now  be  used  to  model  the  system.  The  system 
shown  in  Figure  3.1  is  a  complete  working  model  of  Differential  GPS  constructed  with 
SIMULINK. 

SIMULINK  is  a  non-linear  simulation  package  based  on  MATLAB©.  This  lan¬ 
guage  uses  block  diagrams  and  embedded  MATLAB  routines  to  model  arbitrarily 
complex  systems.  SiMULiNK  allows  simulations  to  be  carried  out  in  continuous  or 
discrete  time.  It  also  has  several  choices  of  integration  routine  —  Euler,  Adams, 
Adams-Gear,  Runge-Kutta,  etc.  Simulation  results  can  be  viewed  as  they  happen 
or  sent  to  the  MATLAB  workspace  for  future  analysis  or  plotting. 

Figure  3.1  depicts  the  inputs  to  the  DGPS  simulation  —  “Aircraft  position”, 
“Satellite  positions”,  and  “Differential  Station  position”  —  which  are  used  to  com¬ 
pute  pseudoranges.  The  receiver  position  and  the  satellite  positions  are  multiplexed 
together  before  they  enter  the  “Aircraft  Pseudorange”  block  (multiplexing  assembles 
an  arbitrary  number  of  separate  vectors  into  one  vector  for  use  by  a  program  block). 
The  “Aircraft  Pseudorange”  block  produces  pseudoranges  complete  with  ephemeris 
error,  ionospheric  delay,  and  tropospheric  delay.  The  receiver  noise,  clock  error,  and 
Selective  Availability  error,  are  subsequently  added  to  the  pseudoranges.  The  remain¬ 
ing  portion  of  the  simulation  removes  some  of  the  errors.  The  “correct  Pseudoranges” 
block  tfikes  receiver  position,  satellite  positions,  and  corrupted  pseudoranges  as  a  mul¬ 
tiplexed  input  vector,  and  removes  the  non-random  tropospheric  and  atmospheric 


39 


Figure  3.1:  DGPS  Model 


delays.  Lastly,  the  atmospherically  corrected  pseudoranges  are  further  corrected  using 
the  “DELTAS”  from  the  differential  station.  Six  differential  pseudoranges  are  outputs 
of  the  model  and  become  inputs  to  a  Kalman  filter. 

The  calculation  of  the  “DELTAS”  is  accomplished  as  follows.  As  in  the  Aircraft 
pseudorange  calculations  above,  the  “Differential  Station  Pseudoranges”  accept  the 
Differential  Station  position  and  the  Satellite  positions  as  a  multiplexed  vector  input. 
The  pseudorange  is  calculated  with  a  routine  identical  to  that  used  in  calculating  the 
Aircraft  Pseudoranges,  but  with  independent  random  variations.  Corrupted  differen¬ 
tial  station  pseudoranges  au’e  outputs  of  the  “Differential  Station  Pseudoranges”  block. 
Clock  error,  receiver  noise,  and  selective  availability  noise  are  added.  The  “correct 
Pseudorainges”  block  is  identical  to  the  one  used  to  adjust  the  aircraft  pseudoranges. 
Having  used  the  Pythagorean  theorem  to  pre-compute  the  exact  ranges  from  the 
differential  station  to  the  satellites,  the  simulation  now  calculates  the  “DELTAS”, 
the  error  in  the  differential  station  pseudorange.  These  values  are  used  to  adjust  the 
aircraft  pseudoranges. 

Next,  each  of  the  major  sections  of  the  model  will  be  discussed  in  more  detail. 

1.  Ephemeris  Model 

The  “Satellite  Position”  block  computes  the  noisy  space  vehicle  positions 
in  tangent  plane  (Local  Geodetic)  coordinates  for  use  by  the  DGPS  model  (see  Fig¬ 
ure  3.2).  The  satellite  positions  enter  from  the  left  in  earth-centered,  earth-fixed 
Cartesian  coordinates  (Conventional  Terrestrial  System).  The  block  “Convert  to 
latitude,  longitude,  elevation”,  a  MATLAB  m-file  ecef2U.m  which  encodes  the  trans¬ 
formation  defined  in  Subsection  2.,  Section  B.  in  Chapter  II.  ,  must  first  be  used  to 
transform  the  satellite  positions  to  the  geodetic  system.  Next,  the  ephemeris  errors, 
resolved  in  the  geodetic  system,  are  added  to  the  space  vehicle  positions.  Finally,  the 
corrupted  positions  are  transformed  to  the  tangent  plane  system  by  the  block  “Con- 
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Figure  3.2:  Ephemeris  Model 

vert  to  tangent  plane”,  a  MATLAB  m-file  U2tanp.m  which  resolves  the  corrupted 
satellite  location  in  the  local  geodetic  or  tangent  plane  coordinate  system. 

2.  Pseudorange  Computation 

The  pseudorange  computation  blocks  “Aircraft  Pseudoranges”  and  “Dif¬ 
ferential  Station  Pseudoranges”  take  space  vehicle  positions,  receiver  (either  aircraft 
of  differential  station)  position  as  applicable,  take-off  time,  atmospheric  pressure  and 
temperature  and  computes  the  pseudoranges  as  shown  in  Figure  3.3.  The  block  “con¬ 
vert  seconds  to  minutes”  merely  converts  the  output  of  the  clock  in  seconds  to  minutes 
and  seconds.  This  value  is  added  to  the  “take-off  time”  in  order  too  maintain  accu¬ 
rate  time  for  use  in  the  ionospheric  delay  calculation.  The  “Ionospheric  delay”  block 
uses  this  time  and  computes  one  delay  to  be  applied  to  all  six  pseudoranges.  Unlike 
the  model  previously  discussed,  a  generic  set  of  coefficients  (q„  and  /3„)  has  been 
assumed  for  this  model  yielding:  A=  20  nanoseconds  and  P=12  hours.  The  scale 
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Figure  3.3:  Pseudorange  Computation  Model 

factor  SF,  is  assumed  to  be  one.  In  order  to  achieve  the  25%  random  error  discussed 
in  Section  1.,  the  output  of  the  “Ionospheric  delay”  block  is  multiplied  by  a  gain  of 
0.125.  The  output  of  this  gain  block  is  multiplied  by  white,  Gaussian  noise  with  unity 
variance,  by  adding  this  value  to  the  nominal  delay,  the  desired  nominal  ±  25%  delay 
is  achieved  95%  of  the  time.  Since  this  value  applies  to  all  six  pseuodranges,  the  one 
value  is  converted  to  a  vector  of  six  with  the  multiplexer. 
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The  block  “Exact  ranges”  uses  the  MATLAB  subroutine  exact. m  to  com¬ 
pute  the  exact  geometric  range  from  the  receiver  to  each  satellite.  It  accepts  the 
receiver  position  and  the  positions  of  the  six  satellites  as  a  multiplexed  input  vector. 
Using  the  Pythagorean  theorem  in  three  dimensions,  this  block  outputs  the  exact 
ranges  to  all  six  space  vehicles. 

The  reminder  of  the  simulation  diagram  calculates  the  tropospheric  delay. 
The  “Satellite  elevations”  block  accepts  the  receiver  position  and  the  satellite  posi¬ 
tions  and  computes  the  elevation  angle  to  each  satellite.  Since  all  quantities  are  in  the 
tangent  plane  Cartesian  coordinate  system  already,  this  is  a  very  simple  operation. 
First,  the  vector  from  the  receiver  to  the  satellite  is  found  by  subtracting  the  receiver 
position  vector  from  the  satellite  position  vector 

^  ^aat  ^rcvr 

y  ”  Vsat  yrcvr  •  (3‘1) 

-  ^  ^  ^  ^rcvr 

Next,  the  vector  is  normalized  to  unit  length 


Finally,  the  elevation  angle  is  given  by 

0  =  arcsinl.  (3.3) 

Note  that  normalizing  the  vector  and  taking  the  arcsin  is  an  equivalent  process. 

The  “Tropospheric  delay”  block  takes  the  output  of  the  “Satellite  eleva¬ 
tions”  block,  the  surface  temperature,  and  the  surface  pressure  as  inputs  to  compute 
the  tropospheric  delay  to  each  individual  satellite.  This  block  executes  the  MATLAB 
function  tropdel.m  which  is  the  tropospheric  delay  algorithm  discussed  earlier  imple¬ 
mented  in  MATLAB.  The  7.5%  random  error  is  added  by  multiplying  the  output  of 
the  “Tropospheric  delay”  block,  the  six  element  vector  of  delays  to  the  satellites,  by 
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0.075  to  form  a  “vector  gain”.  This  output  is  multiplied,  element  by  element,  by 
white  noise  with  unity  variance  to  produce  a  scaled  noise  vector.  The  result  of  this 
operation  is  the  output  of  the  block  “Element  by  element  multiplication”.  Adding 
this  error  vector  to  the  nominal  tropospheric  delay  vector  yields  the  desired  reuidomly 
varying  delay  vector. 

The  outputs  of  each  of  these  processes  —  exact  range,  tropospheric  delay, 
and  ionospheric  delay  —  are  summed  to  form  the  noisy  pseudoranges.  An  identical 
process  occurs  in  finding  the  differential  station  pseudoranges.  However,  the  noise 
added  in  the  differential  station  calculation  is  uncorrelated  with  that  added  in  the 
aircraft  calculation. 

3.  Clock  Model 

The  clocks  of  the  differential  station  and  the  aircraft  receiver  are  both 
modeled  as  previously  shown  in  Figure  2.6.  The  variances  of  the  white  noise  for  the 
aircraft  and  differential  station  clocks  are  identical,  but  the  random  seeds  are  different. 
This  m2ikes  the  white  noise  in  the  two  models  uncorrelated.  The  outputs  of  the  clock 
models  are  multiplexed  into  six  element  vectors  and  multiplied  by  the  speed  of  light 
to  convert  them  to  range  errors. 

4.  Selective  Availability  Model 

Selective  availability  is  modeled  identically  for  each  satellite.  As  shown  in 
Figure  3.4,  the  error  is  white  noise  put  through  a  low  pass  filter.  As  discussed  in 
Section  2.,  the  low  pass  filter  has  a  time  constant  of  three  minutes  (180  seconds).  All 
of  the  white  noise  sources  driving  the  errors  for  each  satellite  are  independent  because 
they  have  different  random  seeds. 
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Figure  3.4:  Selective  Availability  Model 
5.  Pseudorange  Correction 

In  a  manner  similar  to  the  way  the  pseudoranges  are  calculated,  the  cor¬ 
rupted  pseudoranges  are  corrected  with  known  tropospheric  and  ionospheric  models 
as  shown  in  Figure  3.5.  The  only  difference  between  this  part  of  the  model  and 
the  initial  pseudorange  calculation  is  that  this  portion  subtracts  the  nominal  delays 
in  an  attempt  to  remove  the  atmospheric  errors.  Obviously,  this  correction  cannot 
compensate  for  the  random  part  of  the  errors. 

B.  DGPS  MODEL  VERIFICATION 

This  model  of  Differential  GPS  was  verified  by  simulating  it  under  the  following 
circumstances: 

•  receiver  is  stationary  at  the  origin  of  the  tangent  plane  system 

•  differential  station  is  stationary  at  the  origin  of  the  tangent  plane  system 
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Figure  3.5:  Pseudorange  Correction 

•  six  satellites  are  uniformly  distributed  at  elevation  angles  greater  than  ten  de¬ 
grees 

•  the  origin  of  the  tangent  plane  system  is  at  45°  N  latitude  and  45°  E  longitude 
(on  the  earth’s  surface) 

The  values  of  the  differential  station  exact  ranges  are  computed  directly  (with  the 
Pythagorean  theorem)  in  the  model.  Since  the  aircraft  receiver  is  stationary  at  the 
same  location  as  the  differential  station,  any  differences  between  the  aircraft  pseudor¬ 
anges  and  the  exact  ranges  are  errors.  Comparison  of  the  two  outputs  labeled  “psrng 
aircraft”,  “psrng  aircraft  (corrected)”,  and  “psrng”  (for  only  satellite  ^\)  confirms 
the  validity  of  the  model.  That  is,  the  magnitude  of  the  mean  error  (the  difference 
between  the  aforementioned  outputs  and  the  exact  range)  at  each  of  these  stages 
should  decrease  at  each  subsequent  output.  Examination  of  the  “DELTAS”  output 
reveals  the  benefit  of  the  Differential  correction. 
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1.  Verification  with  Clock  Error 

a.  Atmospheric  Correction 

As  Figure  3.6  shows,  the  atmospheric  correction  subroutine  decreased 
the  pseudoranges.  This  figure  depicts  the  difference  between  the  pseudoremge  before 
it  enters  the  correction  block  and  after  it  exits.  Since  the  ionospheric  and  tropospheric 
delays  are  assumed  to  always  be  positive,  the  correction  routine  should  always  reduce 
the  pseudoranges,  resulting  in  a  positive  correction  on  the  plot.  The  increaise  in  the 
atmospheric  correction  with  time  is  due  to  the  ionospheric  error  increasing  because  of 
the  time  of  day.  The  simulation  starts  at  1300  local,  one  hour  before  the  maximum 
ionospheric  delay.  Thus,  the  correction  monotonically  increases  and  would  continue 
to  increase  until  1400  local. 


Figure  3.6:  Atmospheric  correction  for  pseudorange  (with  clock  error) 
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b.  Differential  Corrections 

All  differential  corrections  increased  during  the  simulation  time  (ten 
minutes),  reaching  about  4000  meters  for  some  of  the  satellites.  Figure  3.7  shows 


Figure  3.7:  Differential  Corrections  for  pseudorange  #1  (with  clock  error) 

the  differential  correction  for  pseudorange  #1  monotonically  increiising  during  the 
simulation.  The  reason  for  the  constant  increase  in  the  differential  correction  is  its 
clock  instability.  The  differential  correction  algorithm  can  only  interpret  differences 
between  the  pseudorange  and  exact  range  to  a  given  satellite  as  error.  Unfortunately, 
this  means  that  any  error  in  the  differential  station  clock  results  in  an  error  in  the 
differential  correction.  In  fact,  this  situation  causes  differential  correction  #1  (as  well 
as  the  other  undepicted  corrections)  to  increase  so  rapidly, 
c.  Differentially  Corrected  Pseudoranges 

This  exacerbating  effect  of  the  differential  clock  error  on  the  pseudo¬ 
range  error  is  further  illustrated  in  Figure  3.8.  While  the  pseudorange  error  is  slowly 
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Figure  3.8:  Corrected  and  differentially  adjusted  pseudorange  errors  to 
satellite  #1  (with  clock  error) 

diverging  before  it  is  differentially  corrected,  it  is  obviously  worse  after  the  correc¬ 
tion.  The  reduced  accuracy  of  the  differentially  corrected  range  is  due  to  the  clock 
divergence  discussed  in  the  previous  section. 

2.  Verification  without  Clock  Error 

The  clock  error  of  the  Differential  station  causes  the  differential  correction 
to  diverge,  however,  the  clock  error  is  removed  by  the  Kalman  filter.  Therefore,  to 
better  understand  the  effect  of  the  differential  correction,  the  simulation  should  be 
repeated  without  clock  errors.  This  should  highlight  the  effects  of  selective  avail¬ 
ability,  receiver  noise,  and  the  differential  correction  and  can  be  achieved  by  simply 
disconnecting  the  clock  error  segment  from  the  rest  of  the  model. 

In  this  “clock  error  free”  model,  one  expects  much  better  accuracy  in  the 
differentially  adjusted  pseudoranges.  Selective  availability  errors  should  be  largely 
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canceled  by  the  differential  adjustment.  The  primary  remaining  errors  are  receiver 
noise,  both  on  the  aircraft  and  at  the  differential  station.  The  mean  pseudorange 
errors  should  be  near  zero;  the  standard  deviation  of  the  error  should  be  one  to  seven 
meters  [Ref.  2,  p.  64]. 

a.  Atmospheric  Correction 

As  expected,  the  atmospheric  correction  remained  the  same  (see  Fig¬ 
ure  3.9).  It  turns  out  that  atmospheric  effects  are  totally  independent  of  the  clock 


Figure  3.9;  Atmospheric  correction  for  pseudorange  #1  (without  clock 
error) 

errors.  Thus,  they  are  expected  to  remain  unchanged  in  this  simulation, 
b.  Differential  Corrections 

The  benefit  of  the  differential  correction  can  be  seen  in  Figure  3.10. 
The  high  frequency  variations  in  the  differential  correction  are  due  to  receiver  noise 
and  ephemeris  error.  The  low  frequency  variation  is  clearly  due  to  selective  availabil- 
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Figure  3.10:  Differential  Corrections  for  pseudorange  #1  (without  clock 
error) 

ity.  Selective  Availability,  being  the  only  low  frequency  error  in  this  model,  must  be 
causing  the  low  frequency  variation.  Since  SA  is  common  to  both  the  differential  sta¬ 
tion  amd  the  aircraft,  the  differential  correction  accounts  for  it.  Therefore,  SA  errors 
are  eliminated  from  the  aircraft  pseudorauge  error  by  the  differential  correction, 
c.  Differentially  Corrected  Pseudoranges 

Figure  3.11  clearly  shows  the  benefits  of  Differential  GPS.  The  dif¬ 
ferentially  corrected  errors  are  now  virtually  zero  mean.  That  is,  the  differential 
adjustment  has  removed  the  biases  from  the  pseudorange  errors.  The  standard  devi¬ 
ations  of  the  differentially  adjusted  pseudorange  errors  are  quite  reasonable.  Table  3.1 
summ^lrizes  the  means  (/x)  and  standard  deviations  (<t)  of  the  differentially  corrected 
pseudorange  errors  for  each  satellite  in  meters.  These  figures  show  that  the  mean 
pseudorange  errors  to  each  satellite  are  practically  zero.  Clearly,  the  differential  cor- 
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Figure  3.11:  Corrected  and  differentially  adjusted  pseudorange  errors  to 
satellite  #1  (without  clock  error) 


TABLE  3.1:  MEAN  AND  STANDARD  DEVIATION  OF  DIFFER¬ 
ENTIALLY  CORRECTED  PSEUDORANGE  ERRORS  (WITHOUT 
CLOCK  ERROR) 


1 

Satellite  1 

Satellite  2 

Satellite  3 

Satellite  4 

Satellite  5 

Satellite  6 

im 

-0.0060 

-0.0133 

-0.0051 

-0.0498 

0.0024 

0.0127 

m 

2.765 

1.831 

1.602 

1.653 

2.456 

2.319 

rection  has  resulted  iu  the  error  becoming  zero  mean. 

The  standard  deviations  are  all  on  the  order  of  two  meters.  Using  a 
representative  PDOP  (Position  Dilution  of  Precision)  of  three,  the  DGPS  fixes  from 
this  system  should  be  accurate  to  within  four  meters  (one  rms).  These  values  are 
consistent  with  the  one  to  seven  meter  values  published  for  DGPS  accuracy. 
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C.  INS  SIMULINK  MODEL 


The  INS  SlMULlNK  Model  is  composed  of  the  “strapdown  IMU”  and  the  “INS 
Computations”  blocks  as  shown  in  Figure  3.12.  At  this  top  level,  the  INS  accepts 


Figure  3.12:  INS  Model 

“vdot”  and  “states”  0,  and  from  which  it  ultimately  computes 

“measured  a_u”  (^a).  The  “strapdown  IMU”  block  computes  exact  thrust  acceleration 
and  subsequently  subjects  it,  as  well  as  the  angular  rates  and  Euler  angles,  to  various 
sensor  errors  and  produces  the  measured  outputs.  The  “INS  Computations”  block 
subsequently  uses  the  sensor  outputs  to  compute  inertial  acceleration. 

These  two  blocks  will  be  expleiined  in  more  detail  in  the  remainder  of  this 
section. 

1.  Strapdown  IMU 

The  “strapdown  IMU”  depicted  in  Figure  3.13  is  comprised  of  four  sub¬ 
groups  —  “Compute  thrust  acceleration”,  “Rate  Gyros”,  “Accelerometers”,  and  “In¬ 
clinometers”.  Modeling  of  the  rate  gyros  and  accelerometers  include  all  of  the  error 
sources  —  biais,  noise  floor,  and  cross-axis  sensitivity  —  defined  in  Section  E.  of 
Chapter  II.  Inclinometer  outputs  are  corrupted  by  coupling  with  linear  acceleration. 
Since  Kuechenmeister  developed  the  IMU  model  in  [Ref.  16],  this  portion  of  the 
model  will  not  be  further  discussed  or  verified. 
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Figure  3.13:  Strapdown  IMU  Model 
2.  INS  Computations 

The  “INS  Computations”  block  converts  the  sensor  outputs  from  the  IMU 
into  inertial  acceleration.  The  structure  of  this  block  is  represented  in  Figure  3.14. 
Inertial  body  acceleration  in  the  body  frame  is  computed  by  the  “compute  a_b”  block 
and  provided  to  the  “convert  a_b  to  a_u”  block  for  transformation  to  the  inertial  frame. 
The  “Estimate  lambda”  block  provides  an  optimal  estimate  of  the  Euler  angles  based 
on  its  inputs,  “lambda”  (from  the  inclinometers)  and  “omega_b”  from  the  rate  gyros. 
These  estimated  Euler  angles  are  used  by  the  “convert  a-b  to  a_u”  block  to  compute 
the  rotation  matrices  for  the  transformation. 
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Figure  3.14:  INS  Computations  Model 

Before  transforming  the  acceleration  from  the  body  to  the  inertial  frame, 
the  iiiodel  must  compute  the  Euler  angles.  This  task  is  performed  by  the  “Estimate 
lambda”  block  shown  in  Figure  3.15.  The  top  half  of  this  figure  provides  measured 
Euler  angles  and  calculated  Euler  rates  for  the  bottom  half  of  the  figure,  the  com¬ 
plementary  Kalman  filters,  to  estimate  lambda.  The  block  labeled  “S(lambda)  * 
omega”  computes  Euler  rates  from  Poisson’s  equation  (Equation  2.33)  as  previously 
discussed.  The  “rearrange  into  rate/amgle  pairs”  block  pairs  each  uuler  angle  with 
its  derivative.  Each  pair  then  enters  its  own  filter  —  “Roll  filter”,  “Pitch  filter”,  and 
“Yaw  filter”  —  which  provides  an  estimate  of  each  respective  angle.  The  result  is  fed 
back  to  be  used  in  the  calculation  of  the  Euler  rates  as  shown. 

Finally,  with  “a_b”  and  “Lambda”  computed,  the  “convert  a_b  to  a_u” 
block  in  Figure  3.14  can  calculate  inertial  acceleration,  “a_u”,  through  the  transfor¬ 
mation  defined  in  Equation  2.34. 

D.  INS  MODEL  VERIFICATION 

As  previously  stated,  the  IMU  Model  will  not  be  verified.  However,  the  remain¬ 
der  of  the  INS,  the  INS  Computations  portion,  must  be  verified.  Since  this  portion 
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of  the  INS  is  purely  computational,  no  particular  flight  regime  need  be  tested.  The 
circumstemces  of  the  verification  simulation  are: 

•  Bluebird  aircraft  flying  at  sea  level 

•  elevator  fixed  at  0° 

•  lateral  controls  fixed  at  0° 

•  trim  thrust  applied 


Figure  3.15:  Euler  angles  estimator 


There  are  two  sets  of  data  to  be  compared  in  this  verification.  The  ‘‘INS  Compu¬ 
tations”  block  computes  both  inertial  acceleration  and  Euler  angles.  By  comparing 
these  quantities  with  the  correct  values  calculated  from  the  uncorrupted  kinematic 
quantities,  one  can  evaluate  the  accuracy  of  the  model.  Thus,  the  remainder  of  this 
section  will  consist  of  comparing  exact  and  calculated  values  of  inertial  acceleration 
and  Euler  angles. 

1.  Euler  Angles 

Figure  3.16  depicts  the  difference  between  the  exact  roll  angle  (which  re- 


Figure  3.16:  Actual  and  estimated  roll  angles 

mained  zero)  and  the  roll  angle  calculated  by  the  INS.  At  first  glance,  the  difference 
appears  to  be  substantial.  However,  when  one  realizes  that  the  maximum  difference 
between  the  calculated  and  estimated  roll  angle  is  a  mere  0.42°,  one  recognizes  that 
there  is  really  very  little  error.  This  small  error  is  due  primarily  to  cross-axis  sensi¬ 
tivity  in  the  rate  gyros.  This  phenomenon  is  causing  a  small  difference  between  the 
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measured  and  actual  roll  rate.  Furthermore,  as  the  phugoid  motion  damps  out,  all 
angular  rates  will  tend  to  zero.  Therefore,  the  cross-axis  error  will  tend  to  zero.  This 
is  certainly  an  effective  estimator  in  the  roll  or  X-axis.  The  estimates  of  the  yaw 
angle  behaved  similarly  and  were  also  quite  accurate,  with  the  maximum  error  only 
0.6°.  The  error  is  again  a  result  of  cross-axis  sensitivity. 

Pitch  angles  estimated  by  the  complementary  filter  were  also  quite  accurate 
(see  Figure  3.17).  The  maximum  error  between  the  estimated  pitch  angle  and  the 


Figure  3.17:  Actual  and  estimated  pitch  angles 


actual  pitch  angle  is  only  1°.  Furthermore,  as  can  be  seen  from  the  figure,  the  error  is 
decreasing  toward  zero.  The  minimum  error  is  driven  by  the  noise  in  the  IMU.  Since 
the  angular  estimator  relies  on  the  virtually  noise-free  inclinometer  at  steady  state, 
the  error  in  the  estimate  at  steady  state  would  certainly  be  nearly  zero. 


2.  Linear  Acceleration 


The  inertial  acceleration  calculation  is  quite  accurate,  just  as  the  Euler 
angle  estimation  scheme  is.  Since  the  accelerometers  and  rate  gyros  are  subject 
to  nearly  identical  types  of  errors,  it  is  not  surprising  that  the  acceleration  errors 
would  be  similar  in  nature  to  the  angular  errors.  For  example,  Figure  3.18  shows 
actual  and  calculated  inertial  X  acceleration.  Once  again,  the  values  are  quite  close. 


Figure  3.18:  Actual  and  calculated  X  acceleration 

with  the  maximum  deviation  being  around  0.5  ft/s^.  Moreover,  the  error  appears 
to  be  decreasing  as  the  phugoid  motion  damps  out.  The  error,  however,  will  not 
go  to  zero.  Because  the  Z  accelerometer  will  continue  to  sense  thrust  acceleration 
in  steady,  level  flight,  the  X  measured  acceleration  will  have  a  small  steady  state 
error  due  to  cross-axis  coupling.  Furthermore,  the  “fuzziness” ,  caused  by  noise  in  the 
accelerometer  measurement  of  thrust  acceleration,  clearly  identifies  the  calculated 
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curve.  As  expected,  this  high  frequency  component  remains  constant  over  the  entire 
curve. 

Figure  3.19  more  clearly  exposes  the  effects  of  cross-coupling  and  noise 
in  the  accelerometers.  Although  the  aircraft  never  actually  accelerates  in  the  Y 


time  in  seconds 

Figure  3.19:  Actual  and  calculated  Y  acceleration 

direction,  cross-coupling  with  both  the  X  and  Z  accelerometers  causes  the  general 
downward  swing  of  the  calculated  curve.  While  this  curve  also  appears  to  be  quite 
noisy,  the  standard  deviation  of  the  noise  is  actually  only  a  fraction  of  a  ft/s^.  Like 
the  X  mecisured  acceleration,  the  error  will  tend  to  a  small,  non-zero  value  because 
of  cross-coupling  with  the  Z  accelerometer.  Even  so,  the  calculation  of  the  inertial 
acceleration  in  the  Y  axis  is  satisfactory. 

Finally,  Figure  3.20  shows  the  accuracy  of  the  inertial  acceleration  calcu¬ 
lation  in  the  vertical  or  Z  direction.  This  acceleration  error  does  not  appear  to  be 
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Figure  3.20:  Actual  and  calculated  Z  acceleration 

converging  to  zero  as  do  the  accelerations  in  the  other  axes.  This  lack  of  convergence 
is  made  more  clear  by  Figure  3  21.  This  effect  is  due  to  the  fact  that,  unlike  the  other 
accelerometers,  the  Z  accelerometer  will  not  measure  at  or  even  near  zero  acceleration 
in  steady,  level  flight.  In  fact,  the  one  ‘G’  of  acceleration  that  the  Z  accelerometer 
senses  in  steady,  level  flight  causes  small  errors  in  the  other  two  axes  due  to  cross-axis 
sensitivity. 

E.  NAVIGATION  (LOOKING  AHEAD) 

The  simulations  conducted  in  the  previous  sections  clearly  indicate  which  errors 
dominate  Differential  GPS  and  INS.  Certainly  the  most  troublesome  DGPS  errors  are 
caused  by  clock  differences.  In  fact,  after  only  ten  minutes  of  simulation  time,  the 
clocks  alone  had  driven  the  pseudorange  errors  for  some  of  the  satellites  to  two  miles. 
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Figure  3.21:  Z  acceleration  error 

Simulating  the  system  without  the  clock  errors  showed  that  the  system’s  accuracy 
would  be  vastly  improved  if  the  clock  errors  could  be  eliminated.  Fortunately,  the  INS 
proved  to  be  quite  accurate.  Thus,  no  INS  errors  need  to  considered  in  the  Kalman 
filter  design.  Therefore,  the  Kalman  filter  design  for  navigation  (see  next  chapter) 
must  include  a  clock  model  of  the  difference  between  the  aircraft  receiver  clock  and 
the  differential  station  clock. 
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IV.  DGPS/INS  INTEGRATION 

Having  completely  defined,  encoded,  and  tested  models  of  Differential  GPS  and 
INS,  the  issue  of  integration  can  finally  be  addressed.  As  mentioned  early  in  this 
work,  the  engineering  device  which  accomplishes  the  task  of  fusing  the  outputs  of 
the  INS  (inertial  acceleration)  and  the  DGPS  (six  pseudoranges)  to  produce  inertial 
position  is  the  complementary  Kalman  filter.  A  description  of  the  general  design 
process  for  a  complementary  Kalman  filter  follows. 

A.  GENERAL  LINEAR  KALMAN  FILTER  DESIGN 

Before  detailing  the  Kalman  filter  design  process,  it  is  critical  to  define  the 
framework  in  which  it  operates.  The  mechamism  which  defines  this  process  is  called 
the  synthesis  model.  The  synthesis  model  models  the  process  whose  states  the 
Kalman  filter  is  designed  to  estimate.  In  general,  this  model  is  non-linear.  Since  this 
is  a  linear  design  process,  one  must  first  linearize  the  synthesis  model  to  produce  a 
standard  linear  system  which  can  be  represented  by  a  state  space 

X  =  Ax  +  Biw  -h  B2U 

z  =  Cx  +  Du  +  u,  (4.1) 

where  x  is  the  state  vector,  z  is  the  measured  output  vector,  w  is  the  process  noise 
vector,  u  is  the  input  vector,  v  is  the  sensor  noise  vector,  and  A,  Bi,  B2,  and  C  are 
the  standard  matrices  defining  a  linear,  dynamic  system.  The  feed-forward  matrix. 
D,  is  assumed  to  be  zero.  Complete  explanations  of  all  linear  systems  concepts  can 
be  found  in  any  basic  control  text,  such  as  [Ref.  17,  Ch.  2].  Typically,  w  and  r  are 
specified  as  zero-mean,  white,  Gaussian  noise  vectors  with  covariances  R^,  and 
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respectively. 

Most  often,  a  Kalman  estimator  is  used  in  estimating  states  which  cannot  be 
directly  measured,  although  they  must  be  observable.  For  example,  suppose  the 
measured  output  of  the  system  defined  in  Equations  4.1,  2,  is  a  scalar  which  is  the 
sum  of  all  of  the  system’s  states  (plus  measurement  noise).  A  Kalman  estimator  could 
take  u,  the  system  input  vector,  and  z,  the  measured  system  output,  and  produce  an 
optimal  estimate  of  each  state.  The  estimator  is  constructed  as  follow's 

X  =  At  +  B2U  4-  H{z  —  Ci),  (4.2) 

where  x  is  the  estimator  state  vector,  H  is  the  Kalman  gain,  and  the  remaining 
variables  represent  the  same  quantities  as  in  Equation  4.1.  Rewriting  Equation  4.2 
by  collecting  Fs  yields 

}  =  {A-  HC)x  +  B2U  +  Hz.  (4.3) 

This  form  shows  that  the  stability  of  the  estimator  is  determined  by  the  poles  of  its 
state  matrix,  A  —  HC.  If  this  matrix  has  all  stable  poles,  the  estimates  are  stable. 
But,  w’ill  they  approach  the  actual  states,  x,  that  they  are  supposed  to  estimate'^  By 
defining  the  error,  x,  as  the  difference  between  the  estimator  states  and  the  actual 
states  and  subtracting  Equation  4.2  from  the  first  of  Equations  4.1,  one  finds 

X  =  {A  —  HC)x  +  Biw  +  Hv.  (4.4) 

This  relationship  demonstrates  that  estimator  stability  alone  results  in  the  error  con¬ 
verging  to  zero,  if  there  is  no  process  or  sensor  noise.  The  filter  cannot  be  designed 
to  remove  the  error  caused  by  these  disturbances.  Thus,  the  accuracy  of  the  Kalman 
filter  is  limited  by  the  tw'o  noise  vectors.  To  complete  the  design,  the  Kalman  gain 
must  be  calculated. 

The  Kalman  gain  matrix  is  computed  by  solving  two  equations  in  sequence  — 
tlie  algebraic  Riccatti  equation  and  the  gain  equation.  The  algebraic  Riccatti  equation 
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solution,  Poo,  is  the  steady  state  covariance  matrix  of  the  error  state,  x,  provided  the 
H  is  computed  using  Equation  4.6.  It  is  (see,  for  example  (Ref.  15,  p.  132]) 

APoo  +  PooA^  -  PooC^RrT'CPoo  +  =  0.  (4.5) 

Having  solved  this  equation  for  Poo,  one  can  use  this  result  to  compute  the  Kalman 
gain  with  the  following  equation 

H  =  -PooC'^R-\  (4.6) 

By  inserting  this  gain  into  Equation  4.3,  the  Kalman  filter  is  constructed.  The  re¬ 
maining  issue  to  be  resolved  is  the  specification  of  the  covariances,  Ru,  and  Ri..  It  is 
the  choice  of  these  two  matrices  that  determines  the  characteristics  of  the  filter.  The 
effects  of  varying  these  covariances  are  discussed  in  the  following  paragraph. 

Increasing  the  sensor  noise  covariance,  Rv,  results  in  decreased  gain,  as  one 
can  see  from  Equation  4.6.  Since  the  gain  is  a  measure  of  the  emphasis  given  to 
the  measurement,  z,  this  makes  sense.  That  is,  a  noisier  measurement  should  be 
deemphasized  by  the  filter.  Conversely,  decreasing  the  measurement  noise  covariance 
will  increase  the  gain.  In  the  limiting  case,  as  the  sensor  noise  covariance  approaches 
infinity,  the  filter  will  be  identical  to  the  dynamic  system  w’hose  states  it  is  designed 
to  estimate.  Since  the  Kalman  gain  would  be  zero,  the  measurements,  z,  would  be 
ignored  entirely.  Likewise,  if  the  sensor  noise  covariance  becomes  very  small,  the  gain 
would  be  very  large  and  the  system  would  essentially  invert  the  measured  output  to 
find  the  states. 

The  concept  of  “emphasis”  introduced  in  the  previous  paragraph  is  more  com¬ 
monly  referred  to  as  bandwidth.  Increasing  the  “emphasis”  of  the  sensor  measurement 
by  decreasing  the  sensor  noise  covariance  is  equivalent  to  increasing  the  bandwidth 
of  the  estimator  with  respect  to  that  input.  This  concept  forms  the  basis  for  the 


66 


design  of  a  complementary  Kalman  filter.  Strictly  speaking,  the  process  noise  and 
sensor  noise  covariances  must  be  measured  characteristics  of  the  noise  inputs  to  the 
system.  Realistically,  this  measurement  procedure  is  impractical  and  unnecessary.  By 
using  the  bandwidth  concept,  one  can  tailor  the  estimator  by  selecting  the  covariance 
matrices  to  produce  the  the  desired  frequency  characteristics. 

The  procedure  for  designing  a  complementary,  linear  Kalman  filter  is: 

•  construct  the  synthesis  model 

•  linearize  the  synthesis  model  (if  necessary) 

•  set  both  and  to  the  identity  matrix  of  the  appropriate  dimension 

•  compute  the  Kalman  gain 

•  construct  the  complementary  filter  using  Equation  4.3 

•  determine  the  sensor  bandwidths  using  Bode  analysis 

•  adjust  Ry  to  increase  or  decrease  the  bandwidths  as  required  and  iterate 

This  procedure  was  used  to  develop  the  complementary  linear  Kalman  estimator, 
v.hich  integrates  inertial  acceleration  and  six  DGPS  pseudoranges  to  yield  inertial 
position. 

B.  DESIGN  OF  THE  COMPLEMENTARY,  LINEAR,  KALMAN  FIL¬ 
TER  FOR  DGPS/INS  INTEGRATION 

Before  describing  the  exact  design  process,  the  goal  of  the  process  must  be 
established.  The  final  filter  should  rely  on  the  DGPS  pseudorange  errors  at  low 
frequencies,  below  0.5  Hz  or  about  three  radians  per  second.  At  these  low  frequencies, 
the  filter  should  essentially  invert  the  six  pseudorange  errors  to  find  position.  At 
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higher  frequencies,  the  filter  should  rely  on  the  inertial  accelerations  for  positioning: 
by  integrating  these  accelerations  twice,  the  filter  will  calculate  position.  Based  on 
these  considerations,  the  Bode  plot  of  the  transfer  function  from  pseudorange  error 
as  input  to  corresponding  pseudorange  error  estimate  &s  output  should  look  like  a 
low  pass  filter.  Conversely,  the  frequency  response  from  the  accelerometer  inputs  to 
pseudorange  error  estimates  should  be  similar  to  a  high  pass  filter. 

As  delineated  in  the  previous  section,  the  first  step  in  the  design  process  is  the 
construction  of  the  synthesis  model.  It  is  shown  in  Figure  4.1.  The  synthesis  model 
defines  the  process  whose  states  the  Kalman  filter  will  estimate.  In  this  case,  the 
desired  output  of  the  filter  is  inertial  position  —  P_.x,  P_y,  and  P_z  in  the  synthesis 
model.  Since  each  of  the  scalar  elements  of  the  position  vector  are  states  of  the 
synthesis  model,  they  can  be  estimated  by  the  filter. 

The  basic  construction  of  the  synthesis  model  is  quite  simple.  Each  of  the 
scalar  components  of  the  inertial  acceleration  vector,  a_x,  a-y,  and  a.^,  are  inputs  to 
the  model.  These  accelerations  are  integrated  twice  to  yield  the  inertial  positions, 
P_x.  P.y,  and  P_z.  Three  other  inputs  to  the  upper  part  of  this  diagreim  are  wl, 
w2,  and  w3.  These  are  three  scalar,  independent,  white,  zero-mean,  Gaussian  noise 
inputs.  Integrating  them  before  adding  them  to  the  accelerations  as  sensor  noise 
forces  the  filter  to  compensate  for  an  accelerometer  bias. 

The  bottom  of  the  model  defines  the  unstable  clock  dynamics.  Note  that  this  is 
identical  to  the  previously  defined  clock  model  in  Figure  2.6.  Once  again,  the  purpose 
of  this  portion  of  the  synthesis  model  is  to  alert  the  design  process  to  the  dynamics 
of  the  clock  errors  such  that  it  can  more  effectively  account  for  them.  As  before,  the 
clock  difference  must  be  multiplied  by  the  speed  of  light  to  convert  it  from  a  time 
error  to  a  distance  error. 
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Figure  4.1:  DGPS/INS  Integration  Synthesis  Model 


The  synthesis  model  does  not  contain  two  separate  clock  models  whose  differ¬ 
ence  is  applied  to  the  exact  ranges  to  create  pseudoranges.  The  reason  for  this  is 
quite  simple  —  only  six  independent  integrators  can  be  observable  with  six  measured 
outputs.  Including  two  complete  clock  models  as  well  as  three  independent  process 
noise  inputs  would  result  in  seven  independent  noise  inputs.  Therefore,  the  model 
had  to  be  simplified  to  include  only  five  independent  noise  inputs  in  order  for  all  of 
them  to  be  observable. 

The  “Compute  exact  ranges”  block  uses  the  Pythagorean  theorem  to  compute 
the  ex2u:t  ranges  to  the  six  GPS  satellites.  The  clock  difference  error,  now  converted 
to  distance,  is  added  to  each  exact  range  to  create  six  pseudoranges,  the  measured 
outputs,  2,  of  the  synthesis  model. 

Since  the  use  of  the  Pythagorean  theorem  causes  the  problem  to  be  non-linear, 
the  synthesis  model  must  be  1.  learized  at  a  given  spane  vehicle  constellation  and 
aircraft  position.  The  linearization  process  yields  a  state  space  with  11  states,  five 
process  noise  inputs,  three  control  inputs,  and  six  r.ieiisured  outputs.  P  "  to  the 
complexity  of  this  state  space  model,  the  complel  *  model  is  relegated  to  Appendix 
A. 

The  designer  can  now  determine  the  Kalman  gain.  Unlike  the  general  case,  the 
process  noise  covariance  matrix  is  not  identity  in  this  case.  Since  the  clock  noise 
variances  are  so  small  and  measurable,  they  must  be  specified.  If  the  designer  uses 
unity  variance  instead  of  these  miniscule  values,  the  resulting  filter  will  not  perform 

satisfactorily.  The  process  noise  covariance  matrix  used  in  this  design  is 

■  1  0  0  0  0 

0  10  0  0 

R^=  0  0  I  0  0  (4.7) 

0  0  0  4  X  10-*®  0 

000  0  1.58  X  10-*« 

For  the  first  iteration,  the  sensor  noise  covariance,  Rt,,  is  set  to  a  6  by  6  identity 


70 


matrix  (/  €  /?*).  The  Kalman  gain  is  now  calculated  by  first  solving  the  algebraic 
Riccatti  equation  (Equation  4.5)  and  then  computing  the  gain  (Equation  4.6). 


The  state  space  of  the  filter  is  now  calculated  according  to  Equation  4.3  with 
computed  Kalman  gain,  H.  The  frequency  response  of  the  filter  from  measured 
pseudorange  error,  Ap,  to  the  corresponding  estimated  pseudoranges  error,  Ap,  can 
now  be  evaluated.  In  evaluating  the  bandwidth  of  each  of  the  pseudorange  inputs, 
one  must  consider  the  relative  bandwidth.  Since  none  of  the  zero  frequency  gains 
are  unity,  the  bandwidth  is  three  decibels  less  than  the  steady  state  gain.  The  state 
space  of  the  system  whose  frequency  characteristics  must  be  investigated  is 

i  =  {A-  HC)x-\-  HAp 

Ap  =  Cx.  (4.8) 

After  several  iterations,  the  design  was  finalized  with  the  following  sensor  noise  co- 
variance  matrix 

0.11  0  0  0  0  0 

0  1.1  0  0  0  0 

0  0  0.28  0  0  0 

0  0  0  1.1  0  0 

0  0  0  0  0.11  0 

0  0  0  0  0  0.11 

Bode  plots  from  each  sensor  input  to  its  corresponding  estimate  are  shown  in  Fig¬ 
ures  4.2  through  4.4.  In  all  of  these  plots,  the  ‘‘daish-dot”  line  represents  the  3  dB  less 
than  the  DC  gain.  Thus,  the  intersection  of  the  two  lines  occurs  at  the  bandwidth. 
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Figure  4.2:  Frequency  response  from  pi  to  and  p2  to  p^ 


Figure  4.3:  Frequency  response  from  p3  to  ps  and  p^  to  p^ 
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Figure  4.4:  Frequency  response  from  ps  to  p5  and  pe  to  pt, 


From  the  Bode  plots,  one  can  see  that  all  sensor  bandwidths  are  near  three 
radians  per  second,  as  shown  in  Table  4.1. 

TABLE  4.1:  SENSOR  BANDWIDTHS  FOR  DGPS/INS  KALMAN  FIL¬ 
TER 


Pi 

P2  P3 

P* 

Pb 

pe 

3  dB  bandwidth 

3.10 

3.15 

3.00 

3.15 

3.10 

. 

3.15 

It  is  also  worthwhile  to  note  that,  as  previously  mentioned,  the  frequency  re- 
s[)onse  of  the  transfer  function  from  the  accelerometers  to  the  pseudorange  estimates 
has  negligible  gain  at  low  frequencies.  A  few  of  these  are  shown  in  Figures  4.5 
through  4.7. 
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Figure  4.7:  Frequency  response  from  ^a.  to 

Once  the  the  frequency  characteristics  sought  at  the  outset  of  the  process  are 
aciiieved.  the  filter  can  be  tested  in  simulation. 

C.  FILTER  VERIFICATION 

The  SiMULINK  representation  of  the  simulation  used  to  verify  the  filter  is  de- 
I)icted  in  Figure  4.8.  The  “Linear  filter”  block  accepts  nine  inputs  —  three  com¬ 
ponents  of  inertial  acceleration,  “a.u”,  and  si.x  pseudorange  errors,  “psrng”  - —  and 
provides  estimates  of  three  dimensional  position.  It  is  important  to  note  that  this 
filter,  si'.ice  it  is  linear,  must  be  driven  with  pseudorange  errors,  rather  than  pseudo- 
ranges.  themselves. 

In  order  to  show  that  this  navigation  system  is  functioning,  it  was  simulated 
both  with  and  without  the  accelerometer  inputs.  Since  in  normal,  steady-state, 
trinimed  flight,  the  accelerometer  inputs  are  ignored  (be.,  “washed-out”),  it  is  nec¬ 
essary  to  include  some  high  frequency  dynamics.  By  using  a  five  degree  elevator 
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Figure  4.8:  Integrated  Navigation  System 

half-doublet  (zero  degree  elevator  from  zero  to  ten  seconds,  five  degree  nose-down 
elevator  from  ten  to  20  seconds,  and  zero  degree  elevator  from  20  to  30  seconds),  suf¬ 
ficiently  high  frequency  dynamics  are  excited  enabling  testing  of  the  complementary 
K^dman  filter. 

Due  to  computer  limitations,  the  navigation  model  was  not  be  tested  over  a 
long  period  of  time.  The  computers  being  used  to  perform  the  simulation  did  not 
have  sufficient  speed  or  memory  to  run  the  simulation  for  ten  minutes  of  real  time, 
as  was  done  in  verifying  the  DGPS  Model.  This  shortcoming  prevented  a  detailed 
study  of  the  individual  effects  of  the  INS  and  DGPS  on  the  position  output.  The 
goal  of  this  verification  was  to  show  that  both  sensors  effect  the  output.  This  was 
shown  by  demonstrating  a  difference  between  the  position  errors  with  amd  without 
the  accelerometer  inputs.  One  would  expect  the  position  to  be  less  accurate  with  the 
accelerometer  inputs  since  they  are  so  much  more  corrupted  by  the  various  errors. 
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Both  simulations  were  performed  with  the  following  conditions; 

•  Bluebird  aircraft  flying  at  sea  level 

•  lateral  controls  fixed  at  0° 

•  trim  thrust  applied 

•  r)®  elevator  half-doublet  applied 

1.  Verification  without  Accelerometers 

The  position  estimates  from  the  simulation  without  accelerometers  are 
fairly  accurate.  This  system  navigates  entirely  by  the  differential  GPS.  The  following 
three  plots,  Figures  4.9  through  4.11,  show  the  difference  between  the  estimated  and 
actual  position  in  each  direction  of  the  tangent  plane  system.  The  estimated  position 
does,  in  fact,  track  the  actual  position  quite  accurately  in  all  three  directions.  This 
proves  that  the  Kalman  filter  is  capable  of  effectively  inverting  the  six  pseudoranges 
to  piovide  inertial  tangent  plane  position  without  the  assistance  of  the  INS. 

2.  Verification  with  Accelerometers 

Having  shown  the  performance  of  the  system  with  DGPS  only,  one  can 
compare  the  filter’s  performance  without  accelerometers  to  its  performance  with  ac¬ 
celerometers.  Figures  4.12  through  4.14  show  the  displacement  errors  in  all  three 
directions  of  the  inertial  tangent  plane  system. 

Once  again,  it  is  clear  that  the  position  estimate  tracks  the  actual  position. 
In  fact,  it  is  quite  difficult  to  see  a  great  deal  of  difference  in  the  plots  representing 
navigation  with  and  without  accelerometer  inputs.  However,  reducing  the  data  to 
means  and  standard  deviations  for  each  simulation  reveals  differences  in  accuracies. 
Table  4.2  summarizes  the  difference  between  the  two  simulations  where  is  the 
mean  and  cr  is  the  standard  deviation.  As  forecast,  the  accelerometer  augmented 
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Figure  4.9:  X  displacement  error  without  accelerometers 


TABLE  4.2:  MEAN  AND  STANDARD  DEVIATION  OF  X,  Y,  AND 
Z  ESTIMATED  POSITION  ERRORS  WITH  AND  WITHOUT  AC¬ 
CELEROMETERS 


meters 

X  error 

Y  error 

Z  error 

H  w/o  accel  0.0200 

-0.4481 

0.0607 

<7  w/o  accel  6.323  5.194 

14.84 

//  with  accel 

1.679  0.4093  2.0084 

a  with  accel 

6.594 

5.174 

14.86 

positioning  data  demonstrated  a  slightly  greater  standard  deviation  than  the  pure 
DGPS  positioning.  Furthermore,  it  has  a  significant,  non-zero  mean  error,  while  the 
system  which  did  not  use  the  accelerometers  did  not. 
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V.  CONCLUSIONS 


It  has  been  shown  that  raw  Differential  GPS  outputs  (pseudoranges)  and  raw 
INS  outputs  (inertial  accelerations)  can  be  effectively  blended  to  yield  a  highly  accu¬ 
rate  position  using  a  complementary  Kalman  filter.  Furthermore,  complete  models 
of  INS  and  DGPS  have  been  developed  and  rigorously  verified.  These  high  fidelity 
models  can  be  “transplanted”  to  any  other  work  in  this  area. 

Despite  the  successes  just  mentioned,  there  yet  remains  a  great  deal  of  work  to 
do  before  this  concept  can  be  actually  implemented  on  an  aircraft.  The  additional 
refinements  required  are 

•  discretize  the  entire  system 

•  use  an  extended  Kalman  filter 

•  account  for  changing  satellite  geometries 

•  use  carrier  phase  DGPS 

Each  of  these  elements  will  be  discussed  briefly. 

A.  DISCRETIZE 

All  of  the  modeling  done  in  this  thesis  has  been  done  in  continuous  time,  since 
it  is  really  a  proof  of  concept.  Since  computers  cannot  operate  with  zero  sampling 
time,  this  is  not  realistic.  The  entire  system  must  be  implemented  in  discrete  time. 
For  the  DGPS  model,  this  fact  necessitates  no  change.  Nor  does  it  require  any  change 
in  the  INS  model.  Because  there  is  no  dynamics  in  either  the  DGPS  or  INS  model, 
the  continuous  models  is  identical  to  the  discrete  models. 
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The  dynamic  model  of  the  aircraft  must  also  be  discretized.  This  could  be 
easily  accomplished  by  first  linearizing  the  model  around  a  cruise  condition.  This 
linear  model  could  then  be  e«isily  discretized  by  converting  the  plant  matrix.  .4,  and 
the  input  matrix,  to  their  discrete  time  counterparts.  The  output  matrices,  (' 
and  D  do  not  change. 

Obviously,  changing  to  discrete  time  will  require  the  complementary  filter  to  be 
implemented  discretely  as  well.  This  will  require  the  solution  of  the  discrete,  rather 
than  the  continuous,  Riccatti  equation.  This  equation  is  far  less  computationally 
intensive  and  will  cause  still  greater  rewards  in  the  next  section. 

B.  EXTENDED  KALMAN  FILTER 

The  navigation  system  designed  in  this  work  proved  quite  accurate,  but  it  suffers 
from  some  unrealistic  assumptions.  One  of  these  assumptions  is  stationary  GPS 
satellites.  This  assumption  allowed  the  complementary  Kalman  filter  to  be  designed 
for  a  single,  hypothetical  case  of  space  vehicle  constellation  to  aircraft  geometry.  Since 
the  satellites  are  in  12  hour  orbits,  it  will  be  critical  to  account  for  their  motion  in 
the  actual  system.  Therefore,  an  extended,  complementary,  Kalman  filter  will  be 
required.  The  extended  filter  differs  from  the  regular  filter  in  that  it  is  redesig.ned 
each  time  step.  Rather  than  having  a  fixed  linear  system  that  the  filter  operates  on. 
the  system  is  re-linearized  every  time  step  at  the  estimated  state.  The  new  state 
space  is  used  in  the  discrete  Riccatti  equation  to  find  new  Kalman  gains.  The  new 
state  space  and  Kalman  gain  matrix  is  used  for  only  one  iteration.  Clearly,  with  the 
flight  times  of  more  than  30  minutes,  this  extended  filter  will  prove  necessary. 

C.  ACCOUNT  FOR  CHANGING  SATELLITE  GEOMETRIES 

A  separate  but  significant  problem  resulting  from  the  use  of  stationary  satellites 
is  that  the  effects  of  altering  satellite  to  aircraft  geometry  cannot  be  established. 
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Dilution  of  Precision  is  a  very  troublesome  aspect  of  GPS  and  must  be  considered 
in  a  real  world  navigation  system.  There  are  several  ways  one  could  investigate  the 
effects  of  changing  geometry.  One  way  is  to  actually  model  the  satellites'  orbits  as  a 
junction  of  time.  To  do  this  accurately  could  prove  quite  cumbersome.  The  fact  that 
the  space  vehicles  broadcast  16  different  coefficients  which  describe  their  trajectories 
indicates  the  complexity  of  the  orbits. 

Rather  than  develop  a  model  of  the  entire  constellations  positions  as  a  function 
of  time,  several  geometries  with  stationary  satellites  can  be  investigated.  Since  the 
exact  paths  of  the  satellites  is  of  no  concern,  there  is  no  benefit  to  determining  it 
exactly.  One  can  just  as  easily  verify  that  the  Kalman  filter  is  effectivel}'  accounting  for 
Scitf'llite  geometry  by  merely  considering  several  different  cases  of  satellite  orientation. 
This  method  should  provide  an  equally  sound  test  of  the  entire  scheme  while  greatly 
simplifying  the  process  itself. 

D.  CARRIER  PHASE  GPS 

Even  Differential  GPS  based  on  pseudoranges  turned  out  to  be  insufficiently 
accurate  for  autoland.  If  this  is  the  best  sensor  available,  how  can  this  be  improved? 
I  he  answer  lies  with  Carrier  Phcise  GPS.  This  system  uses  a  phase  difference  of 
arrival  technique,  like  OMEGA,  to  evaluate  position.  Unlike  OMEG.A  with  its  long 
wavelengths,  GPS  signals  have  wavelengths  less  than  a  meter.  In  fact,  static  Carrier 
Phase  GPS  positioning  accuracies  on  the  order  of  one  centimeter  have  been  achieved 
[Ref.  2,  p.  64].  While  the  computational  demands  of  this  may  still  be  too  great  for 
most  computers  light  enough  to  fly,  this  method  may  be  feasible  in  the  near  future. 
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APPENDIX  A:  LINEARIZED  SYNTHESIS 

MODEL 


This  appendix  includes  the  state  space  of  the  model  derived  by  linearizing  the 
non-linear  synthesis  model  depicted  in  Figure  4.1.  The  standard  state  space  matrices 
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Note  that  the  fourth  state  is  clearly  the  clock  difference  while  the  last  three 

states  are  obviously  the  position  states. 
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APPENDIX  B:  MATLAB  FILES 


This  appendix  is  composed  of  all  of  the  MATLAB  m-files  which  are  used  by  the 
SiMULINK  diagrams.  They  are  in  alphabetical  order,  since  no  other  order  suggests 
itself. 

AB-AU.M 

function  au»ab_au(x) 

X  this  fmction  converts  the  noisy  "sensed"  body  acceleration 
X  to  inertial  coordinates  and  adds  gravity  back  in 
X  assign  local  variables 
ab*x(l ;3) ;Eulers=x(4:6) ; 

X  transform  to  inertial  coordinates 
auhat*[l  0  0;0  -1  0;0  0  -l]*ru2b(Eulers) ’*ab; 

X  add  inertial  gravity 
au«(auhat-[0  0  32.174]’); 

ECEF2LL.M 

function  ll=ecef 211(w) 

X  converts  from  earth-centered,  earth-fixed  Cartesian  coordinates 
X  to  latitude,  longitude,  altitude  (geodetic) 
x=w(l) ;y»w(2) ;z=w(3) ; 

X  define  semi-major  md  semi-minor  ellipsoid  aoces 
a«6378 137 ; b»6356000 ; 

X  define  auxiliary  parameters  e  and  f 
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f«(a-b)/a;e»f*(2-f) ; 

X  convert  to  geodetic  latitude,  longitude,  height 
lambda«atan(y/x) ; 

phi*atan(inv (l-e“2) *z/ sqrt (x“2+y“2) ) ; 

X  calculate  N 

N»a/sqrt(l-(e*sin(phi))“2) ; 

X  find  geodetic  height  h 
h»sqrt (x“2+y“2)/cos(phi)-N; 

11= [phi ; lambda ; h] ; 

ECF2TAN.M 


function  tanp=ecf2tan(w, phi, lambda) 

X  converts  from  earth-centered,  earth  fixed  Cartesian  coordinates 
X  to  tangent  plane 

X  given  the  latitude  and  longitude  from  the  workspace 
X  convert  tangent  plane  origin  to  radians 
phi=phi ; lambda=lambda; 

X  transformation  matrix 
T« [-sin (lambda)  cos (lambda)  0 

-sin(phi)*cos(lambda)  -sin(phi)*sin(lambda)  cos(phi) 
cos (phi) *cos (lambda)  cos (phi) ♦sin (lambda)  sin(phi)]; 

X  convert  to  tangent  plane 
tanp=T* (w-112ecef ( [phi , lambda, 0] ) ) ; 
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ELEMULT.M 


function  pro<iuct*elemult(x) 

X  multiplies  15X  of  tropospheric  delay  noise  by  univariant  white  noise 

X  for  randomness  simulation 

for  i»l:6 

new(i)*'x(i)*x(7) ; 

end 

product=new; 

EXACT.M 

fiinction  r  =  exact  (x) 

X  this  function  computes  the  exact  remge  to  six  satellites 
X  given  receiver  position  and  the  satellite  position 
satpos®zeros(6,3) ; 

X  define  receiver  position 
recpos=x(l :3) ; 

X  define  satellite  position  matrix 
for  j®l:6 

satposCj , :)=x(3*j+l ;3*(j+l)) ’ ; 
end 

X  compute  ranges  to  each  satellite 
for  j=l:6 

ra(j)=sqrt((recpos(l)-satpos(j ,l))“2+(recpos(2)-. . . 
satposCj  ,2))''2+(recpos(3)-satpos(j  ,3)) *2)  ; 

end 

r=ra; 
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lODEL.M 


function  io,.del  *  iodel(time) 

%  computes  propaigation  delay  length  for  signal  transmitted 
X  from  GPS  satellite  due  to  troposphere  based  on  time  of  day 
if  time  <  900 
delay=5 ; 

elseif  time  >  1900 
delay*5 ; 
else 

delay*5+20*sin( (time-900) ♦pi/1000) ; 
end 

X 

io.del®delay; 

LAMDQTl.M 

function  lam_dot=lamdotl(x) 

X  computes  Euler  derivatives  given  Euler  angles  and  omega 
phi=x(l) ;theta=x(2) ;psi=x(3) ;p=x(4) ;q=x(5) ;r=x(6)  ; 
lam_dot= [p+q*s in (phi ) *t  an  (theta)  +r*cos (phi ) *tan (theta) 
q*cos (phi) -r*sin (phi) 

(q*sin(phi)+r*cos(phi))/cos(theta)] ; 

LL2ECEF.M 

function  ecef®112ecef (x) 

X  converts  from  geodetic  latitude,  longitude,  elevation  to 
X  earth-centered  earth-fixed  Cartesian  coordinates 
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phi»x(l)  ;laaibda«x(2)  ;h*x(3)  ; 

X  define  semi-major  aind  semi-minor  earth  axes 
a-6378137;b*6356000; 

X  define  auxiliary  quantities  f ,e,  and  N 
f»(a-b)/a;e*f*(2-f) ;N=a/sqrt(l-(e*sin(phi))“2) ; 
X  convert  to  Cartesian 
ecef * [ (N+h) *cos (phi ) *cos (lambda) 
(N+h)*cos(phi)*sin(lambda) 
(N*(l-e“2)+h)*sin(phi)] ; 

LL2TANP.M 


function  tanp=112tanp(x,lat ,long) 

X  convert  from  latitude,  longitude,  altitude  to 
X  tangent  plane  coordinates  given  tangent  plane  origin 

X  convert  to  ecef 
ecef*112ecef (x) ; 

X  convert  to  tangent  plane 
tanp®ecf2t2m(ecef ,lat,long) ; 

REARRANGE.M 

function  new=rearrange(x) 

X  rearreuiges  lambda  2md  omega  into  corresponding  pairs 
new*[x(4)  x(l)  x(5)  x(2)  x(6)  x(3)]'; 
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RU2B.M 


function  Rub  »  ru2b(e) 

X  Euler  angle  transformation  from  {U}  to  {B} 

Tpsi«Ccos(e(3))  sin(e(3))  0;  -8in(e(3))  cos(e(3))  0;  0  0  l]; 
Ttheta*Ccos(e(2))  0  -sin(e(2));  0  1  0;sin(e(2))  0  cos(e(2))]; 
Tphi*tl  0  0;0  cos(e(l))  sin(e(l));0  -sin(e(l))  cos(e(l))3; 
Rub*Tphi*Ttheta*Tphi ; 

SATELEVS.M 


function  elevs»satelevs(x) 

X  this  function  computes  the  elevation  angles  in  radians  of 
X  all  six  satellites 

recpos«x(l:3) * ;satpos(l,  :)=x(4:6) ’ ;satpos(2, :)=x(7:9) ' ; 
satpos(3, :)=x(10: 12) ’ ;satpos(4, :)*x(13: 15) ’ ; 
satpos(5,  :)«x(l''*18)  ’  ;satpos(6,  ;)=x(19:21)  ’ ; 
for  i»l:6 

diff»satpos(i, :)-recpos; 

X  normalize  difference  vector  to  a  unit  vector 
dif f  (3)  =dif f  (3) /sqrt  (dif f  ( 1 )  -2+dif f  (2)  ‘ 2+dif  f  (3)  ‘ 2)  ; 

X  the  angle  is  computed 
satelev(i)®asin(diff (3)) ; 
end 

elevs»satelev ; 
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TAN2ECF.M 

function  ecef*tan2ecf(x, phi .lambda) 

%  converts  from  teingent  plane  coordinates 
X  to  earth-centered,  earth  fixed  Cartesian 
X  given  the  latitude  and  longitude  of  the 
X  origin  of  the  tangent  pleme  system 

X  convert  tangent  plane  origin  latitude  and  longitude  to  radians 
phi»phi=  pi/180 ;  lauBbda=lambda*pi/180 ; 

X  define  transformation  matrix 

T»[-sin(lambda)  -sin(phi) ♦cos (lambda)  cos (phi) ♦cos (lambda) 

cos(launbda)  -sin(phi)^sin(lambda)  cos(phi)^sin(lambda) 

0  cos(phi)  sin(phi)] ; 

X  convert  back  to  degrees  for  112ecef 
phi*phi^ 180/pi ;lambda*lambda^l80/pi ; 

X  convert 

ecef =T^x+112ecef ( [phi , lambda , 0] ) ; 

TANP2LL.M 


function  ll=tanp211(x,lat .long) 

X  converts  from  tangent  plane  to  latitude,  longitude, 
X  and  altitude  given  tangent  pleuie  origin 
X  convert  to  ecef 
ecef»tan2ecf (x.lat.long) ; 

X  convert  to  11 
ll*ecef 211 (ecef) ; 
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TIMING.M 


function  Binute«tiffling(x) 

%  computes  number  of  hours  and  minutes  elapsed  since 
X  clock  start  to  be  added  to  take-off  time  for  use  in 
X  computing  ionospheric  delay 
min»z/60 ; 

X  check  to  see  if  an  hour  has  passed 

if  min  >  60 

hour*round(min/60) ; 

else 

hour»0 ; 

end 

minut e*hour* 100+min ; 

TRQPDEL.M 

function  tropdel»tropdel(x) 

X  used  a  complex  model  to  compute  the  tropospheric  delay  for 

X  all  six  satellites  in  meters 

X  first  calculate  "dry"  portion,  90X 

elev=x(l :6) ;temp=x(7) ;press=x(8) ; 

ae»6378137; 

hdry=148. 98* (temp-4. 12) ; 
for  i=l:6 

idry(i)=(l-(cos(elev(i))/ (l+(l-0.85)*hdry/ae)) “2) “ (-0.5) ; 

sdry ( i ) *2 . 343*pr e ss ♦ ( (temp-4 . 1 2 ) / t emp) * idry ( i ) ; 

end 
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X  now  the  vet  contribution 
hwet* 13000; 
for  i»l:6 

iwet(i)«(l-(cos(elev(i))/(l+(0. 15)*hwet/ae)‘2))“(-0.5) ; 

swet(i)*0.2*iwet(i) ; 

end 

X  total  delay  in  meters  is  the  sum  of  the  two 
tropdel»sdry+swet ; 
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APPENDIX  C:  USER’S  MANUAL 


This  appendix  describes  in  detail  the  steps  required  to  open  and  run  the  software 
which  generated  the  final  results  in  Chapter  IV.  A  basic  familiarity  with  MATLAB 
must  be  2issumed.  Additionally,  it  is  assumed  that  the  user  is  already  logged  on  to  a 
SiMULINK  capable  Unix  work  station. 

Before  entering  the  MATLAB  environment,  one  must  change  the  working  direc¬ 
tory  to  the  one  which  contains  all  of  the  code,  ‘“thesis”,  in  this  case.  The  command 
is 

cd  thesis 

If  the  user  is  remotely  logged  on  to  a  work  station,  he  must  set  the  DISPLAY  en¬ 
vironment  variable  appropriately  in  order  to  display  graphics.  The  command  which 
sets  this  variable  to  intrepid,  a  Sparc  2  work  station  in  the  Avionics  Lab  is 

setenv  DISPLAY  intrepid.aa.nps. navy. mil :0 

Now  it  is  time  to  begin  the  MATLAB  session  by  typing 

inatlab4 

Still  two  more  tasks  must  be  accomplished  within  MATLAB  before  the  simula¬ 
tion  can  be  started.  First,  the  numerous  parameters  which  initialize  the  model  must 
be  loaded  from  a  data  file  with  the  command 

load  bluinit 

To  bring  up  the  SiMULINK  diagram,  type  the  command 

birds 
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It  will  take  approximately  15  seconds  for  the  simulation  diagram  to  appear  on  the 
screen.  Once  it  does,  double  click  on  the  icon  which  looks  like  a  clock  at  the  top  of 
the  diagram.  This  will  allow  you  to  see  the  elapsed  time  as  it  slowly  counts  up  in 
hundredths  of  a  second.  Finally,  a  single  click  on  the  “Simulation”  pull  down  menu 
will  reveal  several  choices.  Choose  the  first  of  them 

Start  “t 

This  particular  simulation  takes  several  hours  to  ruu.  li  outputs  “P”,  the  ac¬ 
tual  aircraft  position,  “v_u”,  the  actual  aircraft  velocity,  “time”,  simulation  time  in 
seconds,  and  “Phat”,  the  estimated  aircraft  position  to  the  workspace.  Once  the 
simulation  is  complete,  these  data  can  be  accessed  like  any  other  MATLAB  data. 
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